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Abstract: A framework of formal physical-chemical kinetics with focus in the mathematical model- 
ling of lumped parameters reaction-diffusion systems is presented. Limit steps for mass transport, ei- 
ther located in the systems boundaries or due to bulk diffusion of certain components inside the sys- 
tems, are employed to give a comparative discussion of the dynamic behavior of several simple mod- 
els of open systems. Two cases of anisocoric (with variable volume) systems whose volume depends 
of the system composition are studied. One case is closely related with Perret and Levey’s biophase as 
a model of a growing population, without scale effects. Considering a biophase as a population of 
growing and dividing spherical droplets, a probabilistic derivation of the linear relation between area 
and volume is given. The other case is an anisocoric generalization of a classical open system model 
due to Ludwig von Bertalanffy, to introduce scale effects in the model. Grounded in slow manifold 
theory, a heuristic generalization of the derivation of Bertalanffy’s growth equation for open physico- 
chemical systems with a complex network of chemical reactions and with scale effect is proposed. For 
growing active droplets with porous matrices, the possibility of a limit step for mass transport due to 
bulk diffusion through the pore space is considered. The stability threshold to droplet division is stud- 
ied using a simplified mathematical model, based on the theory of averaged fields on representative 
volume elements in a two phases medium. Filtration pressure, interfacial tension, viscous-plastic de- 
formation and flow of the polymer matrix, and chemical activity inside the droplet are included. An 
approximate analytical formula is derived for the droplet’s critical radius as a function of reaction 
rates, surface tensions of the droplet interfaces and effective transport parameters of the interiors and 
boundaries of droplets and of the environment. These results are criticized in the light of the progress 
made up to the present time, when a joint work of physicists and biologists in Dresden, Germany, and 
in other places and countries, found several simple mechanisms that could explain how droplets might 
have proliferated, growing, and dividing and, perhaps, evolving into the first cell from an early Earth’s 
primordial soup. 


Key words: Mathematical modelling, dynamical systems, bifurcation theory, slow manifolds, singular 
perturbations, nonequilibrium open systems, formal chemical kinetics, reaction-diffusion systems, aniso- 
coric models, biophase, von Bertalanffy’s and van der Vaart’s mathematical models, active droplets divi- 
sion, prebiotic systems, origin of life. 


(1) Introduction 


In the mathematical approach to open physical-chemical systems far enough from equilibrium, with 
chemical reaction networks and mass transport processes between these systems and their environ- 
ments, three main lines of research can be distinguished: 
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(1)- Mathematical models of the genesis of spatiotemporal dissipative structures through processes of 
successive bifurcations in constant volume (isochoric) systems (“problems of differentiation exclud- 
ing growth’). 


(2)- Mathematical models of variable volume (anisocoric) open physical-chemical systems, perhaps 
including the possibility of a spontaneous division of the system when a critical state is reached, but 
without considering explicitly dissipative structures (“problems of growth excluding differentiation’’). 


(3)- The integration of both points of view, that of growth and that of differentiation, in biomathemati- 
cal research related to the emergence and development of forms and functions in organisms. 


The present paper belongs to the second line of research. The focus is placed in the mathematical 
modelling of the dynamics of nonequilibrium open chemical-physical systems whose volumes depend 
on their compositions, emphasizing possible applications to biology. 

An analytical approach is followed, and closed form analytical formulae are derived when possible 
and their implications are highlighted. 

Insisting on this type of analytical approach may seem a little dated. 

However, to improve a model, it is useful to have an in-depth understanding of the content of the start- 
ing model. An inequality in the space of the model parameters that characterizes the same dynamic 
behavior pattern or an equality that corresponds to a bifurcation point are generally understandable, 
eventually with some computational tools applied to analytical formulas. 

In many cases it is possible to avoid the analytical approach and extract the required information from 
the data obtained in digital simulation runs, using software tools for data analysis, make multidimen- 
sional plots and then examine various projections. But, as Fowkes and Mahoney said, if understanding 
is the aim, then if it is at all possible, extract analytical results (Fowkes and Mahoney, 1994). 


As in a previous presentation by the present author and in a recently published paper ', the lowest level 
of description of the open systems considered in the present research corresponds to a generalized 
formal chemical kinetics with a phenomenological approach to mass transport. In the description of 
droplet’s instabilities, the lowest level of description corresponds to continuum mechanics. So, the 


' The present paper is (in part) an updated and expanded version of an original presentation, given in the 
1978-Latin American School of Physics, combined with an updated, modified, and expanded version of a 
research article published in 2014, in Spanish. 

The presentation “Open physical-chemical systems with composition dependent volume” (July 31, 1978) 
was part of the seminars given in the framework of the 1978-Latin American School of Physics, Maya- 
giiez, Puerto Rico. Some historical and epistemological aspects related with the mathematical description 
of growth in open systems were reviewed. Then, some of the main tools of formal physical-chemical kinet- 
ics necessary to study the dynamics in reaction-diffusion systems were employed to briefly discuss aspects 
of open systems kinetics when the volume is a function of composition and the limit step for mass 
transport is located in the system’s boundary. Some results obtained previously (Suarez-Antola, 1976 and 
1977) were presented by the author and then discussed by the participants. 

The peer-reviewed paper (Suarez-Antola, R. “Bertalanffy’s growth equation as an external approxima- 
tion to the dynamics of Bertalanffy's anisocoric model”, Rev. SCP, 19 (1):41-70, 2014) gives a more com- 
plete framework of formal physical-chemical kinetics. The focus is also in reaction-diffusion systems with 
a limit step for mass transport located in their boundaries. The obtained results were applied to several 
simplified mathematical models of variable volume systems. Based on singular perturbation theory, a 
derivation of a Bertalanffy’s growth equation for open physical-chemical systems with a complex network 
of chemical reactions and with scale effect was given in the paper. (A definition of scale effect can be found 
in Appendix 2 of the present article). 
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lowest level of reference here is mesoscopic nonequilibrium thermodynamics (Kondepudi and Prigo- 
gine, 2002). In certain cases, ad-hoc assumptions will be added to construct simplified models in a 
higher level of description. 


Below the level of mesoscopic nonequilibrium thermodynamics we have the levels of nonequilibrium 
statistical mechanics (Zwanzig, 2001) and quantum dynamics (Dimock, 2011). Both afford a detailed 
theoretical foundation to chemical reaction dynamics (Hendriksen and Hansen, 2012). 

A truly molecular approach to the kind of mathematical modelling of nonequilibrium open chemical- 
physical systems as intended here, is not feasible, even relying in complex numerical algorithms to 
produce atomistic pictures of the kinetics of a chemical process (Unsleber and Reiher, 2020). 

Some of the reasons for this will be considered in the discussion and conclusions section. 

The present distance between the microscopic and the mesoscopic levels of description of open sys- 
tems shows even in the terminology employed. 

For example, in the so-called quantum theory of open systems (Rivas and Huelga, 2012) “open” 
means “not isolated” instead of having the usual meaning of “open” in thermodynamics and formal 
chemical kinetics. 


-In the following, after a brief historical review of some biological research on growth and its mathe- 
matical modelling, a framework of formal physical-chemical kinetics with focus in the mathematical 
modelling of lumped-parameters reaction-diffusion systems, is presented. 

-Limit steps for mass transport, either located in the systems boundaries or due to bulk diffusion of 
certain components inside the systems, are employed to give a comparative discussion of the behavior 
of several simple models of open systems. Bulk diffusion problems are simplified here applying an 
approximate method introduced by Nicholas Rashevsky (Rashevsky, 1960) to allow the derivation of 
closed form analytical formulae. 

-Two cases of systems whose volume depends on the system composition are studied. 

One case is closely related with the biophase (Perret and Levey, 1961 and 1966) as a model of a grow- 
ing population, without scale effects (so system’s surface area is proportional to system volume). 
Considering a biophase as a population of growing and dividing spherical droplets, a probabilistic 
derivation of the linear relation between area and volume is given (Sudrez-Antola, 1977 and 2014). 
The other studied case is an anisocoric generalization of a classical open system model (Bertalanffy, 
1950; Sudrez-Antola, 1977 and 2014) to introduce scale effects in the original model (now system’s 
surface area is proportional to a fractional power of system’s volume). Bulk diffusion effects not in- 
cluded in the above-mentioned references are considered now the biophase model and in what we call 
the “von Bertalanffy’s growing droplet”. 

-A heuristic derivation of a simple growth equation for open physicochemical systems with a complex 
network of chemical reactions and with scale effect is given following the guidelines of the already 
mentioned paper in Spanish (Sudrez-Antola, 2014). 

-For growing active (that is, with chemical reactions inside) droplets with porous matrices under filtra- 
tion pressure, the stability threshold of droplet’s original form is considered using a simplified mathe- 
matical model based on averaged field theory. An approximate analytical formula is derived for the 
droplet’s critical radius as a function of reaction rates, surface tensions of the drop interfaces and 
transport parameters of the drop and its environment. 

-To end, a critical appraisal and final discussion of the above-mentioned results is done in the light of 
recent progress in the research of prebiotic systems and the origin of life. 

-Besides, in Appendix 1, we give a brief introduction to some of the main results about deterministic 
and stochastic dynamical systems than have interest in the framework of the present research. We 
emphasize three methods applied to reduce the number of state variables: Inertia, Centre and Slow 
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manifolds (Suarez-Antola, 2019). 

Some properties of both finite dimension positive nonlinear and positive linear systems, that corre- 
spond to the mathematical models of open physical-chemical systems studied in this paper are briefly 
reviewed. 

In Appendix 2 we give the definition of objects of the same shape as combinations of star-shaped re- 
gions and a brief analysis of the relation between these objects and the scale effect. 

In Appendix 3 we summarize some reciprocity theorems of continuum mechanics that allows the deri- 
vation of an average relative rate of change of a convenient dimension of a viscous body, in terms of 
volume and surface forces. This result is applied to study possible losses of stability of active colloidal 
droplets in section 7. 


(2) Historical aspects 


The growth of many organisms presents a dependence on time such that when the organism's mass or 
volume is plotted in ordinates and time in abscissa, a typical sigmoid curve is obtained: growth is slow 
at first, it accelerates remarkably later , and finally slows down as much as at the beginning and ends 
up approaching asymptotically a maximum volume or mass. 

At the beginning of the last century, based on what was known about chemical kinetics in closed sys- 
tems with constant volume, it was assumed that these characteristics could be explained by a limiting 
step consisting of an autocatalytic reaction. 

Robertson's book summarizes the efforts made in that direction (Robertson, 1923). 

An autocatalytic reaction can be described, from the perspective provided by formal chemical kinetics, 


by means of a logistic differential equation in terms of the concentration x(t) of the reaction product 


(as usual f represents time instants) , a kinetic constant k and a constant c, equal to the sum of the 


initial concentrations of reactant and product: 


—=k-x-(c,—x) (2.1) 


When the initial value x(0) is small enough (x(0) < ah x(t ) tends to c, following a sigmoidal 


Cc 
trend as desired. When x(0) < oi the sigmoidal behavior is particularly evident. 


In 1929, Snell studied a closed-system mathematical model, with an autocatalytic reaction and a vol- 
ume that increased linearly with the amount of product. He concluded that the hypothesis of an auto- 
catalytic master reaction did not explain the observations on the growth of organisms (Snell, 1929). 


In 1934 von Bertlanffy, based on the work of the German physiologist Potter, and in the framework of 
the then still incipient theory of open systems, proposed the following equation to describe the in- 
crease in mass of an organism associated with growth processes (Bertalanffy, 1960): 


tan SoM (2.2) 


The term 7-S resulting from anabolic processes, is proportional to a certain surface S associated 


with the organism whose growth is investigated, while the term «-M resulting from catabolic pro- 
cesses is proportional to the mass M of the organism. For a constant mass density 9, the mass 


M = p-V is proportional to the volume V of the organism under consideration. 


In principle it seems possible to assume that a relationship of this type is verified between the surface 
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and the volume: 
S=pu-V? 0<d<1 (2.3a) 


In the expression V° the volume is assumed to be divided by a unit volume. 


In general: 
é 
S V 
—=4h°'| = (2.3b) 
So " i 
batt M ee ee Pan 
So, considering that V = —, Bertalanffy's growth equation can be rewritten like this, with q@ = 5! 
Pp 
at = a-M?-e-M (2.4) 
t 


In general, it is observed that the parameter 6 can take values between 0 and 1. The (assumed posi- 
tive) coefficients @ (anabolic constant) and « (catabolic constant) are estimated from the adjustment 
of the corresponding solution of equation (4) to the growth curve of the considered organism (or popu- 
lation), determined from experiments. 7 


In a classic work, published in 1950 in the journal Science, Bertalanffy (Bertallanffy, 1950) analyzed 
the differences between the steady state in an open physicochemical system with constant volume 
(isochoric) and the equilibrium state in a closed physicochemical system, relating the properties of 
open systems with the ability to self-regulation and with the (by the so-called) equifinality that organ- 
isms present. In this model, the open system is assumed to have a constant volume V and a constant 
boundary area S . 

Although we will reconsider later the mathematical model for this system in a more general frame- 
work, we summarize the original version of von Bertalanffy here as an historic reference. 

A chemical species A called nutrient, whose concentration in the environment is X can cross a barri- 
er whose permeability for that nutrient is P . Once inside the system it is uniformly distributed with a 


concentration c, and participates in two chemical reactions: an anabolic and reversible one, that pro- 


duces a species B called constituent that cannot cross the boundary of the system, and an irreversible 
catabolic reaction that produces a species C called excreta that can go out to the environment through 


the boundary of the system. The uniform concentration of constituent is c,. The uniform concentra- 
tions of the excreta inside the system and in its environment are c, and Z respectively. The permea- 


bility of the boundary for the excretais R. 


The kinetics constants for the forward A-— B and backward B > A anabolic reaction are k , and 


2 Tf M (0) is the initial mass of the system, the solution of (2.4) for O0<d<1_ is: 


1 
6 2 
M (t) 23) Oe any (0)? -exp| —(1-6)-«-t| When 6 <1 (for example, when 6 = — ) and the 
K 3 
initial mass is small enough, a sigmoidal growth is predicted by the solution (2.5) with an asymptotic mass 
1 


a |r ; 
M_= [| When 6 =1 and the anabolic constant@ is greater than the catabolic constant K, growth is 
K 


exponential: M (t) =M (0)-exp[(a@-x)-1] 
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k., respectively. The kinetic constant of the irreversible catabolic reaction A—>C is k,. 


The mass balance for this open system is given by the following kinetic equations in terms of the con- 
centrations inside the system and in its environment: 


dc, _S 
ay P(X 6) +(-h Gg thy.) +(-k-4) (2.1 a) 
Oe = ky “kG (2.1 b) 
dc, _S 
ay R(e-a)+h-4 (2.1 ¢) 


This simple mathematical model will be generalized in section 4, to include a variable volume de- 
pendent of the composition and to consider bulk diffusion effects inside the system and in the envi- 
ronment. 


If the organisms divide repeatedly, as occurs in a bacterial population when the environment presents 
(and maintains) adequate characteristics, the mass or total volume of the population may increase, as 
time passes, exponentially. 

In 1960 Perret (Perret, 1960) published an extensive verbal, non-mathematical analysis of the growth 
of a bacterial population considered from the perspective provided by the theory of open physico- 
chemical systems capable of expanding (anisocoric), introducing the concept of biophase and linking 
it to the origin of life. 

Subsequently, in 1961 and 1966, Perret and Levey analyzed a mathematical model of a biophase with 
a single chain of first-order, non-catalyzed reactions. In parallel, some works were published on exper- 
imental models of inorganic open systems with catalyzed reactions, the volume of which was varied in 
a controlled way. 

The Perret and Levey biophase can be imagined as a population of homogeneous droplets (spherules), 
whose contents are perfectly mixed, which grow, divide and inside which the chemical reactions that 
make up the mentioned chain take place. The barrier that stops the biophase from its environment is 
permeable to the first metabolite (nutrient) but impermeable to all the others, which are retained inside 
the system. Nutrient flow is assumed to be proportional to the concentration difference between the 
environment and the interior of the biophase and to the separation surface between the biophase and its 
environment. The volume of the biophase is assumed to be proportional to the mass of the last metabo- 
lite in the chain. 

Perret and Levey introduced an ad-hoc relationship of proportionality between the surface and the 
volume of the biophase. When the external nutrient concentration, which is a model parameter, ex- 
ceeds a threshold, the biophase expands asymptotically as an exponential, like a bacterial population in 
the so-called logarithmic growth phase. The exponential coefficient of the exponential expansion in- 
creases with the difference between the nutrient concentration and its threshold value. When the exter- 
nal nutrient concentration is below the threshold, the biophase shrinks and disappears. 


The mathematical model of an open system with constant volume due to Bertalanffy (Bertalanffy, 
1950) and the mathematical model of an open system with an ad-hoc relationship of proportionality 
between the surface and the volume of a bio-phase, due to Perret and Levey (Perret and Levey, 1961), 
were reviewed in a book of thermodynamics and kinetics in biochemistry (Bray and White, 1966). 
Bray and White posed the problem of including scale effects in the mathematical models of open sys- 
tems with variable volumes. (These authors denominate scale effect to the decrease of the surface / 
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volume ratio that occurs when the system volume increases.) 

The scale effect was introduced and studied in a generalized anisocoric version of Bertalanffy’s iso- 
choric model by the present author, in two booklets published in Spanish (Suarez-Antola, 1976 and 
1977) 

The linear relationship between surface and volume, assumed by Perret and Levey in their model, 
requires that the component droplets of the biophase divide asynchronously when reaching a critical 
size. In 1977 (Suarez-Antola, 1977) a detailed probabilistic foundation of the afore mentioned linear 
relationship between surface and volume appeared, based on the strong law of large numbers. 

In an article in Spanish (Suarez-Antola, 2014) this foundation was improved, and some of its conse- 
quences in relation to growth models were examined. This probabilistic foundation of the linear rela- 
tionship between surface and volume is reproduced in the present work. 

Open systems with exponentially variable volume and exponentially variable amounts of reactants but 
constant reactant concentrations (a pseudo steady state) were studied by several authors (Grainger and 
Gaffney,1966) (Grainger and Gaffney, 1968). 

The work of Feinberg, Horn, and others (Feinberg and Horn, 1974) on the formal kinetics in continu- 
ous flow, well stirred and isochoric isothermal reactors allowed to relate the algebraic-topological 
properties of the chemical networks with stability of multiple steady states and necessary conditions 
for sustained oscillations in the concentrations of some chemical species. 

Later, the thermodynamic constraints in nonequilibrium networks of chemical reactions were studied 
for reactions of biological interest in isochoric open systems both with deterministic (Schnakenberg, 
1977) (Feinberg, 1980) (Ederer and Gilles, 2007) (Craciun et al, 2019) and stochastic (Qian, 2016) 
methods. 


A whole cell modelling framework, intended to describe growth and division from the point of view of 
biochemical networks in open system with its volume given as a function of its composition was con- 
structed (Morgan, Surovtsev and Lindhal, 2007). These researchers assumed a constant osmotic pres- 
sure inside a spherical cell, given by the ideal gas law, and a limit step to mass transport located at the 
cellular membrane, considered as an additional phase with constant thickness, besides the phase 
formed by a homogeneous cell interior. The volumes of the two phases change at different rates. As 
consequence, certain oscillations in the components of the cell interior are produced. Although the 
authors introduce a distinction between a phase of cellular growth and a phase of cell division, this is 
done adding further ad-hoc assumptions. 


In the twenties and thirties of the last century the stability of fluid droplets began to be investigated. 

It was carefully studied by experimental means, mainly by physicochemists doing research in the col- 
loidal state (Frey-Wyssling, 1953) (Shaw, 1992). 

Some research was done with active droplets inside which chemical reactions took place. Spontaneous 
drop division was found in certain cases (Fox, 1978) (Evans and Wennerstr6m, 1999) and mathemati- 
cal models intended to explain this behavior were developed and applied prematurely to explain bio- 
logical cell division (Rashevsky, 1960). 

The development of research in cell biology, including advances in the mechanisms of cell division, 
made these pioneering efforts inapplicable. 

However, research on the stability of fluid droplets continued and interest in the possible division 
mechanisms of prebiotic systems reinforced interest in these topics (Whitfield and Hawkins, 2016) 
(Seyboldt and Jiilicher, 2018). 

Besides, there are new technological advances involving interfaces in fluid dynamics which boosted 
basic research in fluid dynamics, including droplet formation through flow instabilities (Gallaire and 
Brun, 2017). 
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To end this partial historical review, let us remind that at near the end of the sixties van der Vaart re- 
considered the classical logistic growth models. He built a new mathematical model of an open and 
constant volume (isochoric) system with an irreversible autocatalytic reaction, leaving the analysis of 
some properties of his model for future work (van der Vaart, 1968). 

However, he says that, depending on the parameter values, the concentration of the product (one of the 
two state variables in this case) could present either a transitory phase with an excess with respect to 
its final stationary value, or damped oscillations as it tends towards the final value. These results sug- 
gest that the old conjecture about the need to introduce an autocatalytic master reaction to explain the 
sigmoid growth curves should be discarded. 

The present author could not find any evidence that van der Vaart published the analysis he had prom- 
ised. So, an ab-initio derivation of some of the results mentioned by van der Vaart was done (Suarez- 
Antola, 2014). This is the main content of the next section, intended as a complement to the historical 
review. 


(3) Results related with the isochoric model of van der Vaart 


In the framework of a critical analysis of autocatalytic growth models, van der Vaart built a model of 
an open system with two chemical species A and B. 

He assumed that both species could cross the boundary of the system, proportionally to the concentra- 
tion difference between the environment and the interior of the system. 

The boundary was assumed to be limit step for mass transport. 

In the interior of the system there is an irreversible autocatalytic reaction: 


A+B52B+C (3.1) 


Instead of writing the kinetic equations as usual for an open system with a constant volume, we will 
write them in a way that can be applied to the variable volume systems of interest to this paper. 


So, let n, represent the number of moles of A and n, the number of moles of B inside the system; J, 


and J, the flow densities of A and B respectively, through the system boundary; g, and g, the 


rates of variation of the amount of A and B respectively, per unit volume, due to the totality of the 
chemical reactions in which these substances are involved inside the system (only one reaction in this 
particular case); S the area of the boundary and V the volume of the system. 

Then we have the following the balance of mass equations: 


dn 
Th =S-J,4V-8, (3.2 a) 
dn 
M2 =$-J,4V-8, (3.2 b) 


The stoichiometric numbers in the irreversible reaction (3.1) for A are v, =0 and v - =1 (asa 


product and as a reactant, respectively) while for B are v Pe =2 and v vs =—1l, soif k is the ki- 


netic constant of the irreversible reaction: 


gi(rm,)= (Cl) (3.3a) 
ga(tm,)=(2-1)-k- he atk tt (3.3b) 
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Representing by P the permeabilities and by X the concentrations in the environment, the flow densi- 
ties are given by the following expressions: 


n 
J =P (x 4) (3.4a) 
Ny 
J,=P,:| X,-—= (3.4b) 
V 
As the volume and the system area both remain constant, we introduce the combined parameter: 
S 
Q=—-P (3.5) 
V 


The concentrations c = as can be taken as new state variables. Then, from equations (3.2) (3.3) (3.4) 


and (3.5) we obtain: 


d 
= = f(c,,c,)= O,-(X,-¢,)-k-¢,-¢, (3.6 a) 
d 
=F lent, )= O, (X, —€, +k +c, +e, (3.6 b) 


When PF. =P, =0 the system is closed and equations (3.6a) and (3.6b) reduce, respectively, to 


dc dc ; 
—!=-k-c,-c, and —*=+k-c,-c,,$80 C,+C, =C, is constant. 
1 2 dt 1 2: 1 2 0 


dt 
7 dc, ae 
From C, +C, =C, and oa =+k-c,-c, an equation equivalent to the logistic equation (2.1) is ob- 
dc 
tained: =e ‘Cy =k-¢,+(Cy-C,) (3.7) 


If the system remains open, the steady state concentrations C,,C, verify, being k>0: 


fi(c,e.) =O, -(X,-¢,)-k-¢,-c, =0 (3.8a) 
fy (c,,c,)=Q, -(X, -¢,)+k-¢,-c, =0 (3.8b) 
Now, we suppose that both permeability parameters are positive. As consequence both combined 


parameters Q, and Q, are also positive. 


From (3.8a) this hyperbolic relationship between the concentrations follows: 


os #4) =2-( 4-1] (3.9a) 


k 1 
From (3.8b) De (15 b) a second hyperbolic relationship is obtained: 
Q y 
2K 
c,=y(¢)= (3.9b) 
C 
k 


However, only non-negative concentration values have physical significance. 


As consequence we to study the functions d(c, ) and y (c,) for non-negative values of the independ- 


ent variable c, and of the dependent variable c, . From (3.9 a) and (3.9 b) it follows that both c, and 
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c, must be positive. Figure 3.1 shows the intersection of the curves c, = o(c,) and c, =y (c,) in 


the point (Z, ; c) that corresponds to the steady state of the system of differential equations (3.6). 


Figure 3.1 


Each function o(c,) and yw (4 ) has two branches. Introducing the restriction c, > 0 in equation (3.9 
a) we obtain a branch of c, = o(c,) that decreases monotonically from +0 approximating to 0 as c, 
grows from 0 to X, after which it becomes negative, so the section corresponding to c, > X, must be 
discarded, leaving 0<c, < X,. 


With the restriction c, > 0 in equation (3.9 b) we obtain a branch of c, =y (eq) that is non negative 


Q 2, 
k 


if and only if O<c,<—. Therefore, only the section with 0 <c, <— = should be retained and the 


k 


rest of the branch discarded. The other branch is always negative, so it must be discarded. 


As the section of Cc, = o(c,) for c, between 0 and X, decreases monotonously (and strictly) from 


+0 to 0, while the section of for c, =y ( A when c, is located between 0 and o , grows monot- 


onously (and strictly) from X., to +00, both non-discarded sections intersect at a single point of 


steady state c,,c, with physical significance. 


Since in this case the Hartman-Grobman theorem (Coddington and Levinson, 1955) (Attle-Jackson, 
1989 and 1991) (Nicolis, 1995) can be applied, the local stability of the point of rest (steady state) of 


the nonlinear system “ =f (c,,¢,) ; ~ =f, (c, Cy) can be studied reducing the study to the local 


stability of the point of rest of the linearized system around that point. 
Let us introduce the new state variables €,=c,—c, and &¢,=c,—C, taking as reference the steady 


state concentrations and define the coefficients of the linearized system in terms of the partial deriva- 
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tives at that point: 


_ AE.) _ AE.) _ f,{C,.€) _ A(E.2) 
aes a fir = fa = foo = 
Oc, Oc, Oc, Oc, 
Then the linearized system becomes: 
d 
On of Bf, (3.10a) 
dt 
ds. 
Hi bit Sons (3.10b) 
Changing to the new variables u,=¢, and u,=f,,:¢,+f\.°¢, the linearized system (3.10) re- 
d d 
duces to: =u, Pe ah 0 1 ab, 
dt dt 
Equivalently, this equation is obtained: 
d* d 
si +2-1- A+ A-u, =0 (3.11) 
dt dt 
1 
Here, by definition:  A=fii-fay—fioe fy (3.12 a) emerscit +f)  (3.12b) 


If A>O and A>O equation (3.11) is the equation for a damped harmonic oscillator and the steady 
state point of the nonlinear system is a locally stable focus. 


From equations (3.8a) and (3.8b) and taking into consideration that 0 <c, < & we obtain: 


fy = Aleut) QO, -k-c, <0 (3.13a) 

C 

yee Weve) 4 a, zip (3.13b) 
Cy 

fy = Fe) 4.5, >0 (3.13c) 
C 

fy = Teo) 9, + 4-6, <0 (3.13 d) 
2 


From (3.13a), (3.13d) y de (3.13b) it follows thata > 0. 
From (3.9a) and (3.9b) this inequality is derived: 


Ong 
d x Ce 2 
— 9(¢,)= fa Q * <—vy(¢)= fu __k 5 (3.14) 
dc, ip k ¢ de, ib Q, _ 
oo 
Considering the inequality — fa < a , the signs of the partial derivatives that appear in the equa- 
12 22 


tions (3.13) and from the equation (3.12a), we have A > 0. 

Then it can be concluded that the steady state of the nonlinear dynamic system (3.8 a) (3.8 b) is a lo- 
cally stable focus: the concentrations approach it by oscillating. 

Van der Vaart's statement about the kinetics of his open system model with an autocatalytic reaction is 
thus supported for the case in which the reaction is irreversible: there is no sigmoid behavior in this 
case, despite the autocatalytic reaction. 
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(4) A generalization of von Bertalanffy open system model and a study of its anisocoric 
versions. 


(4.1) Generalization of the von Bertalanffy’s model to include variable volumes and 
bulk diffusion 


As already said section 2, as a simple example of an open system model capable of self-regulation, 
Bertalanffy (Bertalanffy, 1950) proposed a homogeneous system within which are three chemical spe- 


cies named generically A, B and C. Figure 4.1 shows a draft of this open system (Sudrez-Antola, 
1976). 
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Figure 4.1 


The limiting step for the transport of mass is a barrier located at the border that separates the system 
from its environment. This barrier is permeable to species A and C but it is impervious to species 
B. Once it is inside the system, the species A (called nutrient) reacts reversibly to produce B 
(called constituent) and irreversibly to produce C (called excreta). 

Bertalanffy calls the reaction that connects the nutrient with the constituent anabolic reaction, and the 
reaction that connects the nutrient with the excreta he calls the catabolic reaction. 


He assumes that the three half-reactions are of the first order, with kinetic constants k, (for A— B), 
k, (for BA) y k, (forA>C). 


Since the excreta formation reaction is assumed to be irreversible, this process does not influence the 
dynamics of the nutrient and the constituent. 
As a consequence, the evolution of the state of the system can be described beginning by studying the 


interaction of two variables: n,,n, (the numbers of moles of A and of B respectively). 


If S the area of the boundary and V the volume of the system, a generalized formulation of the dy- 
namic equations of Bertalanffy’s model, that allows to include volume and boundary variations in 
time, is the following: 

dn, dn, 


ShaSS,tV-g, Al-La) “2=V-g, 411b) SAS T4V-g, Allo) 
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Here g, , g, and g, the rates of variation of the amount of A , B and C_ respectively, per unit 


volume, due to the totality of the chemical reactions in which these substances are involved inside the 


system. J, and J, the flow densities of A and C respectively, through the system boundary. 


In case of Bertalanffy’s model, the rates of variation are given by the following equalities: 


n n n 
alms ots)=[ ky Bek, Be} af ky) (4.1.1 d) 
n nN 
Ba(msmy m3) = ky — ky (4.1.1 e) 
nN 
83(m.m, 00, )= hy (4.1.1 f) 


In our present mathematical framework the assumption made by Bertalanffy for the flow densities of 


nutrient and excreta are, respectively J, = P- [x = “| and J, =R- [z = *) In this formulae X 


and Z are the nutrient and excreta concentrations in the environment, while P and & are their re- 
spective boundary permeabilities. 

As the limit step to mass transport is located at the boundary, the concentrations of nutrient, constitu- 
ent and excreta can be assumed to be constant in space inside the system. 

Considering that in principle the volume can vary: 


Be ea (4.1.2a) ajo) (4.1.2b) c,(t) n(t) (4.1.2c) 


V(t) 


However, if bulk diffusion inside the system must be included as a possible limit step to mass 
transport, the geometry of the system must be known to establish the boundary conditions for the 
corresponding partial differential equations in terms of the local concentrations of nutrient, constituent 
and excreta. 

There is a simplified approach introduced by Rashevsky (Rashevsky, 1960) to handle with relative 
ease most diffusion problems (including monoenergetic neutron diffusion in homogenized nuclear 
reactor cores; Suarez-Antola, 2009) making some drastic approximations. 

In our case, let us suppose that is a characteristic length of the system (we can think in an active 


spherical droplet of radius £=7r,) and Dis a representative value of the diffusion coefficient of a 
substance inside the system. If cis a representative value of the substance’s concentration at some 
point between the system centre and system boundary, C, is the substance’s concentration in the sys- 


tem adjacent to the boundary, X is the substance concentration in the environment and P is the per- 
meability of the boundary for the substance, Rashevsky assumes the following approximate equality: 


2-D 
P(X-c,)= F -(cy-e) (4.1.3) 
1 
From (4.1.3): X-q,= -(X-c) (4.1.4) 
£-P 
1+—— 
2:D 
From (4.1.4) we obtain the following approximation for the substance’s flow density: 
P 
fe -(X-c) (4.1.5) 
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measures the relative importance of the resistance to mass 


The positive dimensionless number 


transport of the boundary permeability barrier and the bulk diffusion barrier. When 7 <1 the 


P 
>> 1 the limit step is due to bulk diffusion inside 


limit step is located at the boundary, and when 


the system. This approximate method of handling mass transport can be applied to diffusion in the 
environment (Rashevsky, 1960). 


If the geometric dimensions of the system vary, the characteristic length of the system £ changes and 


P 
the effective permeability (1.2) as well as the substance flow density will be modified by this 
1+ a 


2-D, 


change. 


Let us apply (4.1.5) to the densities of flux of the nutrient and the excreta, assuming that the repre- 
sentative concentrations are given by equations (4.1.2a) and (4.1.2c) and the diffusion coefficients are, 


respectively D, and D, 
Let us suppose that the constituent, besides being unable to cross the system boundary, either doesn ‘t 


diffuse (D, =0) and is uniformly distributed inside the system or diffuses extremely fast that it is 


uniformly distributed also: in any case its concentration is given by (4.1.2b). 
Then from equations (4.1.1) and (4.1.5) we obtain: 


dn, P n 
a (x a) +( ee aes (4.1.6 a) 
ID, 
dn 
et (4.1.6 b) 
dn, _, OR (z tl ky-m (4.1.6) 
dt 14 OR V 
2D, 


The irreversibility of the catabolic reaction A—>C decouples the number of moles of excretan, from 
the evolution equations (4.1.6a) and (4.1.6b) for the numbers of moles of nutrient, and constituent 


n,. However, these last two equations are always coupled and must be studied jointly. 


Equation (4.1.6 a) can be rearranged as follows: 


a ee gO 2 (4.1.7) 
dt l-P V l-P 
1+—— | 2 Sec 
2-D, 2-D, 
Equations (4.1.6) with the addition of a definite relationship between boundary area and system vol- 


ume S (V) and a definite relationship between system volume and system composition V (n, ny, n; ) 
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constitute the generalized version of Bertalanffy’s model of an open system that will be studied here. 


In the original Bertalanffy’s paper (Bertalanffy’s, 1950) the area and the volume of the system are 
constants, and the diffusion processes inside and outside are considered so fast relative to the transport 
through the boundary that they can be neglected . 

e-P C-R 


So both dimensionless numbers and can be neglected relative to 1 and equations (4.1.6) 


reduce to the linear equations of the original model. 


A detailed mathematical analysis of the kinetics corresponding to the original isocoric Bertalanffy’s 
model was carried out (Suarez-Antola, 1976). There is only one globally stable steady state for each 
set of parameters. The concentrations of nutrient, constituent and excreta approximate asymptotically 
to their steady state values from all possible initial conditions. Besides, the concentrations can‘t oscil- 
late during their approximation to their steady values. They increase or decrease monotonically from 
their initial values or present either an overshoot or an undershoot relative to their steady asymptotic 
values. 

In addition, an extension of the model was built in the same publication (Suarez-Antola, 1976), always 
keeping the volume of the system constant. Figure 4.2 shows a draft of the generalized isocoric sys- 
tem. 
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Figure 4.2 
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As shown in the figure, three reaction chains were introduced: a linear chain of reversible reactions 
started in the nutrient, which then branches off into a linear chain of reversible anabolic reactions as- 
sociated with constituents, and into a linear chain of reversible catabolic reactions ending in the excre- 
ta. 

Again, only one steady state existed for each set of parameters values. 

Constructing a suitable Liapunov function and applying the stability theorem of Barbashi-Krasovskii it 
was possible to show the global asymptotically stability of the steady state in the space of states of the 
generalized system shown in Figure 3 (Suarez-Antola, 1976). 


Following the suggestion by Bray and White already mentioned here in section 2 (Historical aspects), 
the model was extended to include the possibility of a composition-dependent volume (Suarez-Antola, 
1977). This approach will be summarized in the next section, but now in the framework of a more 
general model of the open system given by equations (4.1.6). 
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(4.2) Extension of the generalized von Bertalanffy model to include a volume dependent 
on the composition of the system 


Let us assume that the relationship between the surface and the volume is given by equation (2.3a): 


S=pu-V° 0<d<1 (4.2.1) 


Let us further assume the following relation between the system volume and the number of moles of 
the constituent (Suarez-Antola, 1977) : 


V=—-n, (4.2.2) 
Cg 


This imply that the concentration Cg of the constituent remains constant, the same in each point inside 


the system and in each instant of time. 


If the constituent doesn’t diffuse, a constant cg could be compatible with a uniform porous matrix 


formed by the constituent, that allows the diffusion of both, nutrient and excreta, in an aqueous solu- 
tion that fills a fully interconnected pore space through the whole matrix. 


Now, define an effective permeability that combines the effect of the bulk diffusion of the nutrient 
with its permeability due to the boundary barrier: 


Py, =>—————~ (4.2.3 a) 


If the volume varies as a function of the number of moles of constituent, the characteristic length of 
the system ¢ should vary also as an increasing function of this state variable. An approximate relation 
could be the following: 


a Gere 
C(n,)*V3 =| (4.2.3 b) 
Ce 
From (4.2.3 a) and (4.2.3 b): 
P P 
P;(n;) & (4.2.3 c) 
ce 3 
1+ des 1+ ec) 
2-D, ZED Nee 


Substituting (4.2.1) , (4.2.2) and (4.2.3 c) in equation (4.1.6 a) we have: 


O 
an =P, (m){2) [* —Ce tn! —(k, +k,)-n, +k, -n, =f, (n,,n,;X ) (4.2.4a) 
dt = Ny 
Equation (4.1.6 b) remains unchanged: 
Oe hyn, hy = f(s X) (240) 


The only parameter that appears explicitly in the functions f, and f, is the external concentration of 


nutrient X , because, as will be seen below, can be considered a convenient bifurcation parameter (a 
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control parameter as explained in the Appendix 1) to study the changes in the qualitative behaviour 
in the phase space (n, ; Nn, ) ; 
The dynamics associated with these equations has physical sense in the first quadrant of 


the(n,,n, ) plane only, where n, 20 y n, 20. 


2 
Two cases should be distinguished: 0 <6 <1 on the one hand (in particular 6 = s if we think, for 


example, in an active spherical droplet) and 6 =1lon the other. A probabilistic foundation for the 
linear relationship between surface and volume will be developed in section (4.3) below for a popula- 


tion of growing and asynchronous dividing active droplets. 
2 


In Appendix 2, a derivation of the relation S = - V3 for increasing objects of (what is defined there 


as) a constant shape can be found. 


To simplify the analysis, in all parts of this subsection (4.2) we will assume a constant effective per- 
meability P, . This can be a fairly good approximation if the resistance due to bulk diffusion is small 
enough relative to the resistance to mass transport located at the system boundary, or if a constant 
representative value ¢, of the characteristic length of the system / can be introduced to estimate P, , 
if the system is a population of growing and dividing active droplets. 


The dividing droplet volumes are roughly comprised between an initial value v, and an final value 


2-v,, so the characteristic lengths that correspond to the increasing volumes are comprised between 


1 
an initial value (,; *v,* anda final value 0, ~ (1.26)-¢ ;- So, to take an intermediate (, value as a 


L 


representative one to estimate the effective permeability P, = Coe makes sense. 
-P 
1+ 
2-D, 


(4.2.1) A model of a biophase derived from the anisocoric generalization of von Bertalanffy’s 
model of open system. 


P 
When 6 =1, that is, when the scale effect is absent, and P, =7————~_ is an estimation of the 


L .P 
1+— 


effective permeability of the nutrient, the system of differential equations is linear and is a particular 


case of continuous positive linear dynamic system (see Appendix | subsection (9.1.3): 


x 

Baka uet)m al tue Py bm (4.2.5 a) 
@ 

dn 

ne =k,-n,—k,-n, (4.2.5 b) 


There is only one steady state, the origin n, =, = 0. In matrix form: 
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x 
ai" |= —(k, +k, + MPa) Gard * (4.2.5c) 
dt | n, 
k, —k, 


The general solution of the system of equations (4.2.5) can be written as the combination of the two 


1 
linear eigenmodes modes. The first eigen vector k, corresponds to the eigenvalue 2, and the 
A, +k; 
1 
second eigen vector k, corresponds to the eigenvalue A, of the system’s matrix. 
A, +k, 
(1) 1 1 
: n,(t 
In vector form: k 5 = exp(2,(X)-t)-C, ‘|, +exp(A (X)-1)-C, ‘| G2) 
A, +k, A, +k, 
In these expressions C, and C, are constants that are determined by the initial conditions. 
Defining: 
a, =k, +k, +k, + uP, (4.2.7a) 
-P.-k 
ge Xj (4.2.7b) 
Ce 
k k 
Xx. =f Ae | bin, (4.2.8) 
MP, ) k, 


2 
A(x)=-F+ (+) ae (4.2.9a) 

a, a, ; 
A{X)=-> = +a, (4.2.9b) 


From (4.2.9), (4.2.8) and (4.2.7) we obtain: 


2 
a(x)a teehee), [hob eb owl) +H Fok (x-x,) 42.10) 
Ce 


(kK, +4, +k, +p Py) 


lbs 
2 


Cy 


“ky 


A{X)= 


(X-X.) (4.2.10 b) 


[theta 
2 


From formulae (4.2.10) it follows that the first eigenmode is dominant relative to the second. 
While A(x ) is always negative (or with negative real part), A(X ) changes it sign from negative 


(or with negative real part) when X < X. to positive when X > X.,.. 
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As consequence, if X < X. both mole numbers tend to 0, the system consumes itself (is self- 
consuming) and disappears asymptotically when t +o. 


If X > X, the mole numbers grow asymptotically as exponentials with exponent A i( x ) : 


n,(t) *C, -exp| 2,(X)-t] (4.2.11a) 
k 
n,(t)* Tae -C,-exp[ 2,(X)-t] (4.2.11b) 
1 
In vector form this asymptotic behaviour can be written ae * A s exp( P) it Xx) ; t) C1 & 
m\t A,+k 
1 pa 


From (4.2.2) it follows that the asymptotic growth of the volume is also exponential: 


I 
V(t)=—-m(t) ~V,-exp| A(X)-1] (4.2.12) 
@ 


k C 
Here, by definition V, = (4) -—+ depends on the initial conditions through the value of C,. 
TK, ) Ce 


These results show that the nutrient concentration in the environment X is a bifurcation param- 
eter with a critical value X.. 

There are two different patterns of behaviour, one when X < X. (self-consumption) and the other 
when X > X.. (exponentially asymptotic growth). 


The self-consumption scenario corresponds to the global asymptotic stability of the origin, that attracts 
the orbits corresponding to all the possible initial conditions. 

The scenario that presents an unbounded growth corresponds to an unstable origin that repels towards 
infinity the orbits for all initial conditions different to the origin itself. 

The asymptotic exponential growth of the system volume corresponds to the growth of a biophase in 


the sense of Perret and Levey (1961) with exponential growth coefficient A(x ) given by formula 
(4.2.10 a) 


When the nutrient concentration X in the environment is greater than the critical value X. given by 


equation (4.2.8), the exponential growth coefficient increases with X . 
This behaviour is like the one that is observed in the coefficient of exponential growth of a population 
of microorganisms when there is a limiting nutrient. 


d 
However, although for our model of biophase me ) >O is a decreasing function of the envi- 


ronmental concentration of nutrient, A(x ) always increases with X without limit. Figure 4.3 sug- 


gests this behaviour (Suarez-Antola, 1977). 


3 In the original mathematical model of biophase due to Perret and Levey, there are p chemical 
components involved in the linear chain of reactions inside the system. Now, we have p eigenvectors with 
their corresponding eigenvalues, but there is also only one dominant mode that gives the asymptotic ex- 
ponential growth (Perret and Levey, 1966). See Appendix 1, subsection (9.1.3). 
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Figure 4.3 


For a growing population of microorganisms in the so-called logarithmic phase the experimental val- 
ues of A i( x ) in general tend to a maximum when the concentration of the limiting nutrient increases. 
This last behaviour can be recovered if the flow density of nutrient is a suitable nonlinear function of 
X , but if this is the case, the linear model studied in the present section can’t be applied to describe 
the biophase dynamics. 

The nutrient concentration inside the system approximates asymptotically to: 


A(X )+k 
x)= [AG on (4.2. 13a) 
k, 
The constituent concentration is always constant: 
C5= Ce (4.2.13b) 


R 
If Ro = fiat) is an estimation of the effective permeability of the excreta, considering that Z 
eee 


2D, 


is the excreta concentration in the environment, it is possible to show that the excreta concentration 
inside the system approximates asymptotically to + (Sudrez-Antola, 1977): 


A(X )+k -u-R.- 
3,(X)=k,- A(X) tk ; ia sees (4.2.13c) 
A(X)+pu-Re A (X)+k, 


(4.2.2) Scale effects in the anisocoric generalization of von Bertalanffy’s model of open system 


2 
When 6 = 3 , that is, when the system presents a scale effect (Appendix 2), and PR, = P is constant 


because bulk diffusion can be neglected, equation (4.2.4 a) reduces to: 


1 \3 2, 
i y.e[ 1] {x cot ns (k, k,)-n, } ky +n, = f,(m,1;X ) (4.2.14a) 


2 


Equation (4.2.4 b) remains unchanged: 


4 In this reference the diffusion inside the system was neglected, so the formula that appears there for the 


asymptotic concentration of the excretais: AX) =k, (45 )+ = ; ees Tare 
, 1 +H : 1 +k, 


20 


August 2020 


dn, 

dt 
The existence and uniqueness theorem (Coddington and Levinson, 1955) (Markus, 1976) (Attle- 
Jackson, 1989) applies to the system of differential equations (4.2.14 a) and (14.2.14 b) inside the first 
quadrant (n, > 0 and n, > 0), which is an open set G, of the phase plane. 


=k, -n, —k,-n, = f,(n,n,;X) (4.2.14b) 


The steady states (rest points) (7, : ii, ) of the system (4.2.14) are solutions of the equations: 


2 
3 . i 2 
we{t I ioe ae [Eb ).n mS (4.2.15 a) 
Ce k, 
k,n, =k, -n, (4.2.15 b) 
In (4.2.15 a), by definition: 
k 
X,=— Cy (4.2.16) 
k, 


The origin 7, =0 and 7, =Ois always a solution of equations (4.2.15), has physical meaning, but it 
doesn’t belong to Gy. 

If X < X,the origin is the only solution, and as we will see, in this case the system self- 
consumes: from any initial condition (n,(0),n, (0))eG, the state (n, (t),n, (t)) tends to 


(0, 0) and the system asymptotically disappears. 


However, if X > X, there is a second solution of the steady state (fixed points) equations, know be- 


longing to G,: 


2 3 
2-(oe Pk {2} -@-%) (4.2.17 a) 
Cg 
_k 


hain, (4.2.17 b) 
k, 
Now, we will see that for any initial condition (n, (0) Nn, (0)) €G,, different to the rest point (7, ; 71, ) 


given by (4.2.17), the state of the system (n, (t),n, (t)) tends to this non trivial rest point (77,71, ) 


when time increases. 


So, when X > X, and under suitable initial conditions (7, (0) <N,, Nn, (0) <n,) the system pre- 


sents a significant but bounded growth process. 


Then X,, given by equation (4.2.16), is a critical concentration of nutrient in the environment, 


like the critical concentration X , corresponding to the biophase (equation (4.2.8). 


A digital simulation of the dynamics, using the Runge-Ktitta algorithm, provides evidence that con- 
firms these conclusions about the qualitative behavior of the system (Sudrez-Antola, 1977). 
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When the initial volume of the system and the initial amount of nutrient inside are small with respect 
3 
; =~ 1 UP Ek “ es, he ; 
to their stationary values V =— eee (X —X ) and m,=—-n, respectively (being 
Ce k, k, k, 
N, = Cg -V ), the results of the numerical calculation suggest that in general a sigmoid increase in vol- 


ume over time. 


The following table shows a result of numerical computation of system’s volume V (t ) as a function 
of time f. 


The time increment is At=1. Following the differences V (1) -V (t - 1) =AV (t), the sigmoidal 


pattern in volume increase appears: AV (t) increases slowly at first, then accelerates and approximat- 


ing to the steady state it increases less and less. 
The digital simulation was done in1976 (Suarez-Antola, 1977 


t V(t) V(t) -V(t-1) 
0 0,100 
1 0,160. 0,06 
2 0,241 0,08 
3 0,343 0,190 
4 0,467 0,12 
5 0,613 0,14 
6 0,771 0,17 
7 0,965 0,18 
8 1,166 0,20 
9 1,382 0,21 
10 1,610 0,22 
20. 4,03 2,82 
30 5,86 1,83 
40 6,93 1,07 
50 7,47 0,54 
60 7,74 0,27 
70 7,87 9,13 
80 7,93 0,06 
90 7,96 9,03 
100 7,98 0,02 


The following values of the parameters were chosen: k, = 2,k, =k, =P, =M=Cg =1,X =1.5 


The initial values are n(0) =0,0751 and n,(0) =0,1 


As we said in the introduction, our focus here is in the study of qualitative, geometric patterns and 
their changes when a control parameter is modified, deriving analytical formulae for critical values of 
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control parameters and if possible, for the dynamic coefficients that characterize systems behavior as a 
function of time. 


To study the geometric patterns of the orbits, in both cases X > X,and X < X,, we begin by intro- 
ducing the change of variables x=k,-n,—k,-n, and y=n, . 
In the new variables x and y , let us define the new functions f ( y) and g ( y) : 


1 


3 


f(y) =h+h+kht+ uP & (4.2.18 a) 
y3 
A 
g(y)=| ky -k- I y (4.2.18 b) 
y3 
-P-k 
Ae Xa) (4.2.19) 
Ca? 
In the new state variables, the differential equations (4.2.14) transform into: 
d 
—y=a\(x, = 7 (4.2.20 a) 
Bae) 
d 
Oley) = Fy) x8 (9) (4.2.20 b) 


The domain where the theorem of existence and unicity of solutions can now be applied is the open 
half-plane y>Q. This half-plane includes the open set G = {(x, y) ER’ /y>0;x>-k,- y} that 
corresponds to G, in the original positive variables n,,n, . 

When X > X, then from (4.2.19) we see that A>QO . Then the system (4.2.20) has only one rest 
point (0, y.) €G with: 


3 
re k,-k 2 
ks k, +k, +¢ 
2K; ° Ce 
We define: 
y 1 a) 3 
H(y)=], 8(0)-dn=She key ee (4.2.22) 


Figure 4.4. The solid line curve corresponds to A >0 . When A < Othe corresponding curves are the 


specular images reflected in the H (y) axis. 
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1 
Then, multiplying (4.2.20 b) by x and defining V(x, y) = rl +H ( yi we obtain, following a 
trajectory from any initial point (x(0) ,y ( 0)) belonging to the half plane y > 0: 


d 
“V(x y)=-f (9) (4.2.23) 


Considering the definition (4.2.18 a), in (4.2.23) f (y) >k,+k,+k,>0 for any y>0, and as con- 


sequence £V(xy)=-F(»)-* <0 foranyx,and y>0. 


1 1 > 
The function V ( y) = 3 xP + ae -K,° ay = 5 -A-y? is well defined in the whole R’ for any value 


5 2: 


of A, being (-y)3 =-y?. All the level curves V(x, y) =c, are closed. 
They encircle either the point (0, y. ) when A > 0 or the point (0, = y.) when A<0O. 
The level curves form a nested family (filling the whole plane R’) in the sense that if Co > C,, then 


the level curve V (es y) = C,,, is located in the interior of the set bounded by V ies y) =C, >: 


Figures 4.5 (a) and 4.5 (b) show a sketch of these nested curves. 


Figure 4.5(a) A>O Figure 4.5 (b) A<0O 


(A)-Let us first consider what happens in the case X > X,. 


To do this we apply the invariant set theorem (Luenberger, 1979) that appears in Appendix 1, subsec- 
tion (9.1.2) of the present work. 


I 1 3 us 
When X > X,, then A>O and the function V(x, y) aca v5 -ky+y? Se has one sin- 
gle minimum (precisely at the rest point (0, y.)) so (according to subsection (9.1.2) of Appendix 1) 


when A > 0 the function V es y) is a Liapunov function for the system of ordinary differential equa- 


tions (4.2.14) relative to {(x. y) ER /y> o} and corresponding to the fixed point (0, y. ) : 
d 2 te 
Let us introduce the set S = (x, y) EG/ ae (x, y) =F (y) -x° =O? that appears in the invari- 
t 


ant set theorem. S$ = {(0, y) ER*/y> o} corresponds to the positive half-y axis. 
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The only invariant set included in S$ is the rest point (0, y.)s so from the invariant set theorem we 


conclude that all the trajectories beginning in any bounded region Q, = {(x, y) EG/ V(x, y) < s} 
tend to this rest point as time increases. 


As the nested level curves fill R’ (Figure 4.5 (a)), each point es y) €G_ belongs also to one and 


I I Shea, 
only one level curve Vie )H5 > toe ye A ey w TPES, (x, y)EQ,. So, 


given any initial point in G, there is a Q, such that the initial point belongs to it. As consequence the 
trajectory that begins in this point tends to the fixed point (0, y.) : 

From x=k,-n,—k,-n, and y =n, jointly with (4.2.21) and (4.2.17 a) we see that y. =n, and the 
rest point (0, y.) correspond to the rest point (7i,,77, ). As both systems are equivalent in the sense of 


having orbits with the same pattern of qualitative properties ( Appendix 1, subsection (9.1.2)), we con- 
clude that for any initial condition (n, (0),n, (0)) € G, either the state of the system (1, (t),n, (t)) 
tends to the non-trivial rest point (7, ; i, ) when time increases, or it remains there indefinitely. 


. . . “Cc . . * k . 
In short, if the external nutrient concentration verifies the inequality X > X* =—-cg, there is a sta- 
1 


tionary state inside the first quadrant, a point of rest of the dynamic system defined inG,, towards 


which all orbits that start at points of G, tend to. 


(B)-Now, let us consider the other case, when X < X,. 


1 1 3 2 
If X < X,,A<0O and Viay)as to ke ky te Aly? . This function has a single min- 


imum at the point (0, 0). In fact, the origin (0,0) of the x, y phase plane corresponds to the only 
solution with physical meaning of the steady state equations (4.2.15), but it doesn’t belong to the open 
set G={(x,y)eR*/y>0;x>-k,-yt. 


As the nested level curves fill R’ in also this case (Figure 4.5 (b)), each point belonging to G belongs 


I 1 3 2 
also to one and only one level curve of the function V es y) = 3 x + ae -K,° ye + 5 . [| -y?, 


d 
Also, in this case a (x, y) = -f(y) “x” is strictly negative if x #0 and zero when x =0. 
t 


Applying the theorems of regularization and of Vinograd, (that appear in Appendix 1, subsection 


(9.1.2)), the origin (0, 0) becomes a fixed point of a new dynamics that has the same geometric 


pattern of orbits as the original one. From equations (4.2.20) and the transformations given in Ap- 
pendix 1, subsection (9.1.2), the transformed equations are: 


de 1+ Ja? (x, y) +b? (x,y) 1+d?((x,y),{(xy) eR? /y <0}) (ey). 224%) 
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dx 1 d°((x,y).{(xy)eR?/ y <0}) 


dt (1+ ,fa? (xy) +b? (x,y) H(i neeiss) (4.2.24 b) 


However, the invariant set theorem cannot be directly applied because the now fixed point (0,0) be- 


longs to the boundary of the domain G = {(x. y) eR’ /y>0;-k,-y< x} and furthermore, all the 
points of the closed half plane y <0 are fixed points of the new dynamics. 


So, when X < X,,, anew approach to stability seems necessary. 


A weak version of Poicaré-Bendixson theorem for dynamical systems in R’ (Beltrami, 1997) can be 
applied to overcome the mentioned difficulties. To begin with, let us consider the family of closed 


sets, one for everym > 0: C,, ={(x, y)€R?/-k,-ySxSm-ysy >0}U{(0,0)} 


Figure 4.6 shows a sketch of this set and an initial point ( x(0), y (0)) of an orbit, located in the inte- 


rior of C,,. 


| o(x(0), y(O)) 


; x=my 


ent, 


neces, 


DAM “nas sds SO 


an 
| 


My x=-k2.y 
Figure 4.6 


For every possible initial point (x(0), y(0)) eG= {(x, y) ER’ /y>0;-k,-y< x} there is a 
m> 0. such that (x(0), y(0)) E (C,, ), , being (C,, ), the interior of C,, . See Figure 4.6. 


For this same initial point, there is a c, such that if U, = {(x, y) ER*/V (3% y) < c,} (is a closed and 


bounded set) then (x(0), y(0)) E (U, ), . So (x(0), y(0)) ec, qu,. 


The set C,, AU, is closed and bounded, and as consequence is compact. As the only fixed point of 
the dynamics (4.2.24) is (0,0) and belongs to the boundary of C,, 0U,, the weak form of Poicaré- 


Bendixson theorem allows us to conclude that the orbit that begins in (x(0) ,y (0)) tends to (0,0) as 


time increases. 
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Due to the geometric equivalence between the orbit patterns of (4.2.24) and the orbit patterns of 
(4.2.17), beginning from any initial state (n, (0),n, (0)) the open physical-chemical system self- 


consumes and asymptotically disappears. 


Now, let us study in more detail the geometric patterns of the orbits when X > X, and when 
X < X,,, taking advantage of the fact that it is a dynamic system in the plane. 


d. dx, 


IX. 
In an autonomous dynamic system as =f, (a X55 c) mre =f, (ais X55 c) (where c represents a 
t t 


parameter or parameters) time is eliminated and an ordinary differential equation is obtained for 


de, fa (Mx) 
ele 0% of 
F(x, a) c)# dx, FAX: 50) 


The solutions of this last equation (or its reciprocal for f, (x, X55 c) #0) are functions 
xX, = (x; x, (0), c) (or their inverses x, = 9(x;; x, (0), c)) each of which determines an orbit in 


the phase plane. 


Its qualitative properties can be found by applying classical methods in the phase plane (Coddington 
and Levinson, 1955) (Roxin, 1972) (Hirsch, Smale and Devaney, 2004) (Jordan and Smith, 2007). 

The results of the digital simulation confirm the results of the qualitative analysis of the orbits and 
provide information on the trajectories themselves, that is, on how the state variables depend on time. 
The quantitative details of this time dependence are, of course, outside the scope of qualitative meth- 
ods for the phase plane. 


Applying the above-mentioned classical methods, the qualitative characteristics of the orbits in the 


phase plane n,,n, were determined without resorting to numerical calculation. 


In Figure 4.7 the orbits of the equivalent dynamical system (after regularization and application of 


Vinograd’s theorem) are sketched in the half plane n, > 0, for the two situations of interest: X > X, 


and X < X, respectively. 
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Figure 4.7 (a) X > X, Figure 4.7(b) X < X, 


In Figure 4.7 (a) the orbits approaching the origin in the first quadrant of the phase plane change 
their direction and approach to a separatrix curve that connects the origin to the steady state located 
inside the first quadrant. All of them tend asymptotically to this stationary state. 

The diagram below in this figure shows an outline of the mentioned separatrix curve and sections of 
some of the neighboring orbits. This separatrix is both an unstable manifold of the origin and a stable 
manifold of the steady state point inside the first quadrant (Attle-Jackson, 1989) (Roxin, 1972). 


k 
The digital simulation results suggest that when € =— is small with respect to 1, the phase of ap- 


2 
proximation of the orbits to the separatrix occurs in short time scales with respect to the time scales in 
which these same orbits evolve when they are located next to the separatrix curve. 


dn, 


k 

Furthermore, the separatrix curve is located very close to the isocline n, = —>-n, on which We =0. 
t 
1 

In Figure 4.7 (b) all orbits that go through points in the first quadrant approach the origin and tend, 
to this now unique rest point, in the same direction. 
The diagram that appears below in this same figure shows a separatrix, which is approached by the 
other orbits represented there. 


The digital simulation results suggest that, in general also in this case, when € =— is small with 
2 


respect to 1, two-time scales can be distinguished: a short one in the section of those orbits that ap- 
proach to the separatrix, and a long one in the section that tends to the origin in the vicinity of the sep- 
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aratrix. In this case also, the separatrix is located close to the isocline n, = —*-n, in the vicinity of the 
1 


origin. 

When the scale effect is not present, the pattern of the corresponding asymptotic growth dynamics is 
simpler, as we saw in the subsection (4.2.1) “A model of a biophase derived from the anisocoric gen- 
eralization of von Bertalanffy’s model of open system”. 

However, the problem of given an adequate support to the hypothesis of a linear relationship between 
surface and volume in a biophase remains unsolved. This is the subject of the next section. 


(4.3) A probabilistic foundation for the linear relationship between surface and volume 
in a biophase 


The scale effect disappears in a biophase, because in this case 6 =1 so the exchange surface is always 
proportional to the volume of the system. 

To see why this is so, consider a biophase consisting of a set of growing subsystems. 

It is assumed that when one of these subsystems (which can be called spherules) reaches a critical size, 
it becomes instable and fragments into two daughter spherules. The latter begin to grow, until they 
reach a critical size and divide in turn. 

Based on the available experimental results obtained working on populations of proteinoids and other 
analogous colloidal systems (Fox, 1978) (Schrum, Zhu and Szostak, 2010), it is concluded that the 
duration of the division phase of the spherule can be neglected with respect to the duration of its stage 
of growth. So, for the purposes of this paper, division will be considered instantaneous. 


Figure 4.8 shows a representative draft of the volume increase of a microorganism between two cellu- 


lar divisions (mitosis) separated by an interval of time of duration 1, (Rose, 1976). 


A similar pattern of volume change can be observed in the case of an active growing droplet. 
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Figure 4.8 The initial volume is U; and the final volume is U, 


If the biophase originates from the growth and division of a single spherule, a non-negligible correla- 
tion between the times from birth to division can be expected for the first three or four generations, but 
then that correlation can be neglected and the behavior of each spherule, as regards the duration of its 
growth phase, becomes independent of that of the others. 


Then, the random variables $, and ¥, representing the boundary area and the volume of the k-th 


spherule, respectively, can be considered independent in probability, except at a very early stage of 
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biophase evolution. 
Suppose that in an instant of time ¢ there are N spherules. The total area of the biophase at that time 
. WN UN 
is S= a and its volume is V = > Pi ; 
k=l k=l 


The quotient between the area and the volume of the biophase is a new random variable: 


(4.3.1) 


Each of the variables $, and ¥, take values in an interval bounded by a positive lower end (s, for the 
area and v, for the volume of each spherule) and an upper end (s, for the area and v, for the vol- 
ume of each spherule). 

Both s, and v, can be considered as equally distributed with finite variance, being their mean values 
s= E|s| and v=E [7] respectively. 


Then the strong law of large numbers (Welsh, 1970) allows us to conclude that, when N >, 
1 a 
NW -> Vv, converges almost certainly to 


the average volume v and // converges almost certain to: 


N 
Das converges almost certainly to the average area S , 
k=l 


(4.3.2) 


<1 | 41 


H= 


As the scale effect occurs during the growth of each individually taken spherule, it is expected that, at 
each moment of the life cycle of that spherule, between its area § and its volume Y a non-linear rela- 
tionship of the following type is established: 


$x Q(¥)=n,-¥ (4.3.3) 


So, after taking averages: 


Therefore: M= = F My 
v 


When two random variables are related by a in general nonlinear function like § = Av), a Taylor's 
development of ) = V +6,Y around the average V = E [7] gives as a result: 


+ g0(%) -OV+ : aay) (60) +... 


a sees 2 dv 


Averaging S the following expression is obtained: 5 =E [s] = o(v ) + 


LPO) pf (5s) ]+-- 


It can be seen that the approximation S ~ Av ) , which is often done in theory of physical measure- 
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ment errors, requires that the sum of the products of the derivatives of order greater than or equal to 


the second, of the function ov ), with the corresponding central moments of ¥ be neglectable regard- 


ing § > pv). 


However, considering the inequalities s; <8, <s,, v,; SV, Sv, (kK =1,2,...,N ) and the definition of 
ft (equation (4.3.2)), both @ and its almost certain limit 2 can be bounded in this way: 
S. 


ee ee Sy 
sie pe (4.3.4a) ee pet (4.3.4b) 
Ve Vy. Ve Vv; 


If v,= 2-v, and the formula (4.3.3) is applied to connect the area with volume at the same instant for 


the same spherule, from (4.3.4b) the lower and upper bounds that were sought are finally obtained: 


Uo tes 


1— 


1.25 
I 


(4.3.5) 


i i 


With the deduction of the almost certain convergence of the random quotient = to a number 4 


S| NA 


and with the lower and upper bounds given, for this limit, in the inequalities (4.3.5), one of the objec- 
tives of the present work has been reached. 


Next, we study the possibility of approaching the asymptotic growth dynamics of the generalized 
model of Bertalanffy, both with and without (biophase) scale effects, from another point of view. 
This will be the subject of the next section. 


(4.4) A singular perturbations approach to the dynamics: the volume and mass growth 
equations 
As we said while describing Figures 4.7 (a) and 4.7 (b), for the model of anisocoric system with 


2 k 
O= 3 the digital simulation results suggest that when € = — is small with respect to 1, the approx- 
2 


imation of the orbits to the separatrix occurs in short time scales with respect to the time scales in 
which these same orbits evolve when they are located next to the separatrix curve. 


This happens both when X > X, and when X < X,, that is when the separatrix ends in a fixed point 


in the interior of the first quadrant and when it ends in the origin of the phase plane. 
2 
To study the kinetics of the model for 6 = . using digital simulation procedures, it is necessary to 


assign numerical values to seven parameters: the three kinetic constants, the permeability of the barrier 
at the surface of the system, the coefficient that relates the area of the surface to the system volume, 
constituent concentration and the external nutrient concentration. 

The results of the digital simulation runs are interesting and rich in information, but they lack the nec- 
essary generality for the purposes of this work. 
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The purpose of this section is to deduce a differential equation whose solution allows to describe the 
variation of the system volume as a function of time, so that the parameters of this growth equation are 
expressed as a function of the model parameters. 


2 
This will be done here not only when 6 = 3 but more generally for any 6 such that O< 6 <1, and 


for any relationship between nutrient flow density and concentrations in the environment and within 
the system. 

The basis for all this is in the singular perturbation method and slow manifold theory summarized in 
Appendix 1, sub section (9.3.3). 


To begin, let write the equations (4.1.6 a) and (4.1.6b) of the general anisocoric model again: 


Bh = 8B X— BL) (-hym tym) +h) (44.1 a) 
dt V ; 

dn 

Ot = ky, ky ty (4.4.1 b) 


In (4.4 a) appears the effective permeability, introduced in section (4.2), equations (4.2.3), where the 
meaning of the different parameters that appear in the formula are explained: 


P 


ial FCEa 


2-D, 


In the present case the volume of the system could change, during growth, an order of numerical mag- 


(4.4.2) 


nitude or more without dividing, so there is no possibility of substituting (n,) by a representative 


value {, as can be done in the case of a population of dividing spherules. 
In subsection (8.2) of the section of discussion and final conclusions, the consequences of a variable 
effective permeability will be studied. 


But in this section, we take a constant representative value of PR, = P to simplify the analysis. 


It is best to start by replacing the flux density in Eq. (4.4.1 a) with a more general (possibly nonlinear) 
expression that could present the kind of saturation effects that are often observed in microbial physi- 
ology (Rose, 1976): 


n 
h=t(x.%) (4.4.4) 
The following evolution equations are obtained in the phase plane: 
Bis. 6 (x, 2) 4 hm thm) +(-ksm) (4.4.5 a) 
dt V 
Oe = ym ky (4.4.5 b) 
t 
has 1 5 ; 
Considering that V =——-n, and S=y-V° equation (4.4.5 a) can be rewritten: 
Ce 
Sua wvef x, cot + (=k, «1, + k,n) +(—k,-m,) (4.4.5 a bis) 
t ny 
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Introducing the dimensionless variables €& =—~- y &,=—*+, taking into account that 
n n 
i 2 


* 


So ok k . . 
k,n, =k,-n, , X,= ran and V= I Tek &,, equations (4.4.5 a bis) y (4.4.5 b) can be re- 


1 Ce Ce 
written this way: 
dl ye. —ae 4.4.6 
Bo Eat (x, xB} Ch 6-6 )eOb8) (4.4.64) 
oor (E = 4.4.6 b 
wha i 4) (4.4.6 b) 


Defining a new time scale, without dimensions, tf, =k,-tand introducing the parameter ¢ = =, 


2 


equations (4.4.6) are reduced to the following: 


tts (x, x-BheL( k,-(&, + a-( k,-&) (4.4.7 a) 
DB pe 
mas é) (4.4.7 b) 


According to (4.4.7 b), at the beginning ¢, approaches €,, from its initial value ¢, (0), the faster the 
smaller itis €. 
Inspection of equation (4.4.7 a) suggests that €, evolves from its initial value €, (0) more slowly than 
&, when this last variable relaxes from ¢, (0) until it approaches €, . 
Despite the fact that, in general, ¢, (t) it does not tend to é, (t), it gets closer and stays closer, more 
closer the more the smaller is €. 
Therefore, for ¢ small enough with respect to 1, in (4.4.7 b) the approximation €, ~ €, can be made, 
except perhaps in a short initial interval of time (short when compared to the external time scale 
T, = Z ), called the boundary layer or internal scale. 

3 
The theorems that apply in this and other similar cases can be found in the references cited in Appen- 
dix 1, sub section (9.3.3). 


k k 
If —+- (é = €,)can be neglected (which happens if —- is at most of the order of unity) that approxi- 


k, k, 
mation can also be made in (4.4.7 a). So, outside the boundary layer, equation (4.4.7 a) can be approx- 
dé 1 5 1 
imated by: —22. @ ——___y-V°. f (X, X,)+—-(-k,- 
dt, on f ( ) k, ( 3 Gs) 
Returning to the original time we obtain: 
d 1 
Ooo ww L p.v? £ (X, X.)+(-ky-&) (4.4.8) 
dt n, 
or fs Be iis. k, 1 Peas 
Considering that V == = —-- , and — = = , multiplying both members of (4.4.8) 
Ce. - Ce Hele. hotter ae 
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* 


n 
by — finally results in a volume variation equation of the Bertalanffy’s growth type: 
Ce 


dV _ me f(X, X.) 


or = V°-k,-V=a,-V°-k,-V (4.4.9) 


* 


The parameters a, (anabolic) and x, (catabolic) are given, as well-defined functions of the parame- 


ters of the anisocoric open system model, by the following formulas: 


= wie E) (4.4.10 a) ky =k, (4.4.10 b) 
When the flow density of nutrient is linearly related with the difference between the external and the 
S UP, (Xx —X .) 


internal concentrations, Q, . In the nonlinear case, we assume that the flow den- 


xX 


sity of nutrient J, = f (x. X,) is zero when X = X,, negative if X < X, and positive if X > X,, 


increasing with the external concentration of nutrient. 


As consequence, if X > X, the system volume always tends to a non-zero limiting value: 


He. 
V.= (se) (4.4.11) 


4 

If the initial volume (considered in the external time scale, eventually after a fast relaxation in the in- 
ner time scale) is small enough, the system can grow following a sigmoidal pattern. 

If X < X, the volume always tends to zero and the system is self-consuming. 

Therefore, we have rediscovered in this simplified model (4.4.9) the role of external nutrient concen- 


tration X as acontrol parameter. We found the same critical value X, of that parameter and the same 


role in distinguishing two possible dynamic patterns: growth or disappearance of the system. 


The constituent concentrationc,, by hypothesis, remains constant. Nutrient and excreta concentra- 
tions, after the initial relaxation stage on the inner time scale, remain close to and tend toward their 
steady state values. 

So, on the external time scale, the density o of the system can be considered constant. 

Therefore, its mass is proportional to its volume with an invariant coefficient of proportionality on 
the external scale. 

If the parameters a@ = fo ‘ad, and K =k, are introduced, equation (4.4.9) is equivalent to the mass 


growth equation: 


OMY pag hse (4.4.12) 
dt 
But (4.4.12) is an equation identical in appearance to equation (2.4). 
So, by a singular perturbation method, we obtained Bertalanffy's growth equation (2.4) as an approxi- 
mation in the external time scale to the kinetic equations of the anisocoric generalization of the von 


Bertalanffy open system model. Tis is called an external approximation to the original dynamics. 


But there is an important difference: in the growth equation (2.4) the anabolic constant and the cata- 
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bolic constant are empirical parameters, while in the equation (4.4.12) they are given by well-defined 
analytical formulas in terms of the more basic parameters of the model of the physical-chemical open 


1 HOT AG X,) 
xX 


To generalize this approach by singular perturbations and slow manifolds to more complex open reac- 


system: a= (4.4.13 a) K=k, (4.4.13 b) 


* 


tion-diffusion systems whose volumes depend on their composition, we need to extend the framework 
of analysis. 
This will be done in the following section. 


(5) Formal kinetic models for open homogeneous reaction-diffusion systems with com- 
position-dependent volume. 


To make the study of growth (understood as an increase in mass and volume as time passes) inde- 
pendent from the study of differentiation (understood as the genesis and evolution of spatial patterns), 
in this work only homogeneous reaction-diffusion systems are considered. 

The limiting step for mass transport processes is at the boundary of the system or if this does not oc- 
cur, the effects due to diffusion within the system are treated by the approximate Rashevsky’s method 
and a lumped parameters mathematical model is obtained again. 

Furthermore, it is assumed that the systems are isothermal and that electrical phenomena can be ne- 
glected. 

These assumptions will be critically considered in the discussion section at the end of this article. 


The homogeneity hypothesis considerably simplifies the analysis but leaves out of the scope of this 
analysis the questions related to the genesis of possible mechanical instabilities that lead to a sponta- 
neous partition of the system. As consequence, the necessary relationship between the surface and the 
volume of the system must be entered independently. 


The assumption that system’s volume depends on the composition, allows a connection to be made 
with contemporary research, both experimental and theoretical, on that twilight and still poorly under- 
stood region, which separates inert systems from living systems: prebiotic systems. 

Three subsystems interact in the latter: 

(1) A network of chemical reactions, precursor of metabolism. 

(2) A macromolecular element, a precursor to the genome, that carries information: it catalyzes its 
own production and that of the structural elements of the system - the container - regulating the net- 
work of chemical reactions. 

(3) A container, which delimits the prebiotic system of its environment. Able to grow, it regulates the 
transport of mass to and from the environment, ensuring necessary conditions for the operation of the 
chemical network and, where appropriate, the replication of the information-carrying subsystem. 


Having mentioned some of the limitations of the approach, proposed in the present section, to analyze 
prebiotic systems, in what follows we consider an open homogeneous system, consisting of 17 chemi- 


cal species, representing by ”, the number of moles of the i-th species found in the inside the system. 

It is introduced: the area S$ of the boundary that separates the system from its environment, the net 
density of flux (influx less efflux) J; of the i-th species across the boundary, the volume V of the 
system and the rate g, of variation of the quantity of substance i-th, per unit volume, due to the totali- 


ty of the chemical reactions in which this substance is involved inside the system. 
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The balance equation for the number of moles n, can be written like this: 


dn, (dn, dn, : 
rs pee salpeet =S-J,+V-g, i=1,2,...,m (5.1) 
dt dt ), dt “5 
ransporte reacciton 
Since the system is assumed to be homogeneous, the local concentration c, coincides with the aver- 
ee n; 
age concentration — : c,=— i=1,2,...,m (5.2) 


Each elemental chemical hemi-reaction directly connects a complex formed by the reacting species 
together with their corresponding stoichiometric numbers, with a complex formed by the product spe- 
cies together with their corresponding stoichiometric numbers. 

In an elementary reversible chemical reaction two opposite elemental hemi-reactions directly connect 
two complexes to each other. In an irreversible reaction a single hemi-reaction directly connects two 
complexes. 


If r hemi-reactions occur in the system, their hemi-reaction rates are represented by: 
v, =O Gpesge) Prager (5.3) 


It is assumed that the speed of a hemi-reaction depends directly on the concentrations (or activities) of 
the chemical species that form the corresponding complex of reactants, according to the law of mass 
action. 


For this reason, UV, generally depends on a subset of the / species that make up the system. 


To link the algebraic structure of the reaction network with chemical kinetics, the structure is de- 
scribed using a directed graph. 
The p vertices of the graph represent each of the different combinations (complexes) of chemical 


species that appear either as reactants or as products. Each of the Yr directed arcs connect a vertex of 
reactants with a vertex of products when the corresponding hemi-reaction can occur. 
Two matrices and a column vector are introduced: 


-A stoichiometric coefficient matrix A = lv, | of mxp , where v,, is the stoichiometric coeffi- 


cient of the i-th species in the j-th complex. In this approach, stoichiometric coefficients are consid- 
ered non-negative. 


-An incidence matrix E = le ‘ ol of pxr, where e, , takes the value 1 if the  -th arc hits the vertex 


j and is directed towards it (the jth species is produced by the hemi-reaction), and e, , takes the val- 


ue -1 if the o-th arc hits the vertex j but is moves away from this vertex (the jth species is con- 


sumed by the hemi-reaction), and 0 otherwise. 


-A vector P of rx1 is introduced also, whose elements are the weights U 5 (C,,€p5+5C,) of the arcs 


(the po -th component of this vector, U,(c,,C,,...,c,, ) is the expression of the speed of the po -th hemi- 
p p p\C1 22 a p pP p 


reaction). 
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Taking all this into account, the speed of variation g,of the amount of the i-th substance, per unit 


volume, due to the totality of the chemical reactions in which the substance is involved inside the sys- 
tem , can be expressed like this: 


P r 
8 (5Cyl = Se Wess, De (Gtssese,.) (5.4) 


j=l p=l 
Now introducing a vector g of mx1whose components are the overall reaction rates of each of the 
chemical species, equation (5.4) can be rewritten in compact vector form: 


g=AEP (5.5) 


Starting with (5.5), we introduce the vector J of mx1, whose components are the flow densities 
across the system boundary and the vector n of mx1, whose components are the numbers of moles. 
Then the scalar evolution equations in the molar number space can be rewritten as follows, also in 
compact vector form: 
d _ = a = = 
7 =5-3+V-g=S-3+V-(AEP) (5.6) 
t 


p.R 


If the positive numbers v,"" and We 7 represent the stoichiometric numbers of the i-th species when it 


is part, respectively, of the complex of the reactants and when it is part of the complex of the products 
of the hemi-reaction considered, and taking into account that when a species does not appear in a 


reaction complex associated with a given hemi-reaction, its stoichiometric number there is null, equa- 
tion (5.4) can be rewritten as follows: 


GNC Cee) = y’ VP ev (C,ep.00Cy) Hw BORO 77) (5.7) 


p=l 
However, equations (5.4) or their compact version (5.5) are better suited to study, from a mathematical 
point of view, the relationship between open physicochemical systems and the algebraic structure of 
the underlying network of chemical reactions. 


Finally, to apply the mathematical model of formal kinetics to systems capable of expansion, it is nec- 
essary to relate the area of the border to the volume of the system, and the latter to the composition 
(71, ty ye..5My )- 

In the next section, in all cases it will be assumed that the formulae (2.3 a) or (2.3 b) that appears in 
the historical review in section 2, is applied to connect area with volume. 

In section (7) "Active droplet growth and stability" this restriction will be removed 


(6) Heuristic derivation of an equation analogous to the von Bertalanffy growth equa- 
tion 


In section 4 the following equation was obtained (from an equation for volume variation) to describe 
the mass variation from balance between an anabolic term and a catabolic term: 


ah: a-M°-K-M 
dt 
| - | we Peuael BES _ ue f(X, X*) 
The anabolic parameter @ is in turn given by the formulas: @=p " -Q@, and Q, = "ye ; 


where J, =f (x 5X ‘) is the nutrient flow density through the permeability barrier. The flow is can- 
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celed when the external concentration of nutrient X becomes equal to the critical value X,. The flow 


is positive and increasing when X grows from X,. 


Loree , a \I-6 
Under these conditions, if 0 <0 <1 the mass increases approaching a maximum value M, = (<) 


If 6 =1 the mass undergoes exponential growth with a constant 1 =a-kK. 


When X is less than the critical value X,, the flow is negative, and then the terms anabolic and cata- 


bolic are negative: for 0< 6 <1 the system consumes itself and disappears. 
In sum, the external nutrient concentration operates as a control parameter in such a way that the dy- 
namic system presents a bifurcation when that concentration crosses a critical value. 


The deduction of an equation for the growth of volume in the framework of the Bertalanffy anisocoric 
model was made based on the existence of two scales of time, of different orders of magnitude, an 
internal scale and an external scale. 

The rapid relaxation of concentrations, on the internal scale, towards the proximity of certain station- 
ary values made possible a description of the dynamics on the external scale using a single summariz- 
ing variable, the volume of the system. 

Relaxing concentrations to the vicinity of certain stationary values produces in the molar number 
space (state space) a rapid approximation (i.e., on the internal time scale) of the state of the system to a 
manifold known as the slow manifold (because in its adjacencies the state evolves on the external time 
scale). 


The question arises under which conditions this approach could be generalized for a homogeneous 
system with any network of chemical reactions (including autocatalytic reactions). 

In the latter case, the temporal evolution of the number of moles of the component species of the sys- 
tem can be described compactly using the nonlinear vector equation (5.6) that appears in the previous 
section: 


LE pe eer (6.1) 
dt 


Introducing a column vector M whose components are equal to the molar volumes M ; of each of the 


m chemical species that make up the system, the total mass of the latter can be expressed in matrix 
form as follows: 
M(t)=Mii(t) (6.2) 
dM(t) 


From (6.1) and (6.2) we obtain: erm ‘MI (t)+V-M7% 
t 


The law of conservation of mass implies that the mass of the system can only be modified through the 


exchange of matter with the environment, so that M’ g =0 and then: 


dM(t —7= 

aM(t) ie S:-M'J(t) (6.3) 
dt 

The component species of the system can be divided into two groups: the group of constituents I’. 


(which cannot enter or leave the system) and the group of permeants I’, (which can cross the border 


that separates your environment system). 
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The flux density of each permeant species depends on the internal and external concentrations of the 
species considered and possibly on the concentrations of some other species. The latter occurs when 
there are cross-effects in the permeability barrier or modifications in the permeabilities associated with 
variations in internal or external concentrations, or diffusive coupling inside the system. 


The group of permeants is subdivided in turn into the group of nutrients I, and the group of excreta 
| ee 
Nutrients normally enter to the system while excreta normally leave it. 


Taking into account this last partition, the speed of variation of the mass of the system can be rewritten 
as follows, restricting the sums to nutrients and excreta: 


dM 
AMO) 2 YiM,.J,+5->M,.J; (6.4) 
d t icD'y ic, 
When the system grows, in this last equation the nutrient mass flows are all positive and the excreta 
mass flows are all negative. 


The rate of change of the total mass of excreta M,, is given by the equation: 


dM, (t) 
S-YM,.J,+V- YM ,.8, (6.5) 


ieD; ic, 


Now we assume again the existence of two scales of time scales of different order of magnitude, an 
internal scale and an external scale or growth scale. 

Suppose further that during a time interval whose extension is small with respect to the time scale 
associated with growth, the concentrations within the system approach certain stationary values (but 
do not necessarily tend to them, but can oscillate in their environment) and that from that moment the 
increase in the mass of excreta in relation to the increase in the total mass of the system can be ne- 
glected. 

Then in equation (6.5) the approximation can be made: 


S-Y°M,.J,xV->iM..g, 
ieD’, iel', 
From this last equation and from (6.4) a balance equation is obtained between a term associated with 
anabolic processes, proportional to the area of the boundary between the system and its environment, 
and a term associated with catabolic processes, proportional to the volume of the system: 


dM\t — — 
O53 D iy (6.6) 
dt icDy iel', 
Since P= > M, F ie Sif, -c , once the concentrations c, vary close to certain stationary values 
i=l =l 


in , the density can be considered constant during the growth process. 


M 
Then, from V =— and from (6.6), an equation is deduced that corresponds to the original formula- 


ye) 


tion of the von Bertalanffy growth equation (equation (2.2) that appears in the historical review): 


= SM (6.7) 


By definition: 
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n= > M,.J, (6.8 a) 
icDy 
> M..8, 
— ile (6.8 b) 
ve) 


Adding the relationship S = -V° between the area of the boundary and the volume of the system, 


from this last relationship and from the relationship between the mass and the volume, the equation for 
the variation of the mass (2.4) is finally deduced: 


CLE ey, (6.9) 
dt 
q-f 7 (6.10) 
p 


In the subsection (8.2.5) of section 8 “Discussion and final conclusions”, we consider why we can 


» M,.8; 
expect that 77 = > M ;J, and K= a __ are practically constant parameters during system’s 
ily 
growth, although significant fluctuations can be expected in the individual values of flow densities and 
reaction rates. 


However, a basic flaw in all this analysis is the ad-hoc introduction of a relationship between area and 
volume, unavoidable when working with homogeneous systems. 

This flaw will be reconsidered in the discussion section and final conclusions, and it is avoided in the 
next section where the relationship between surface and volume naturally results from the mathemati- 
cal model. 


(7) Active droplet growth and stability 


Now we consider an active spherical droplet of radius % . The droplet material consists in a loose but 


interlinked polymeric matrix with a fully interconnected por space filled by a multicomponent solution 
in water as solvent. 

The other components may diffuse, be advected by a volume flow, and react chemically inside the 
droplet or at its boundary. 

The boundary of the droplet could be a membrane or a simpler interface. 

This system will be homogenized using averaging field theory (Whitaker, 1999) (Quintard, 2015) giv- 
ing a spatial superposition of the matrix with the space of pores. > 

In the next section 8, in the discussion and final conclusions, other and better models of active droplet 
division will be briefly considered. 


5 There are two main approaches to homogenization of polyphasic media: one employs the so-called for- 
mal multiscale asymptotic methods, based in perturbation theory applied to periodic polyphasic media, 
and the other employs volume averaging over a suitable defined representative volume element (a com- 
parison of these two different approaches is given in Davit et al, 2013). The volume averaging method can 
be applied to situations in which an assumption of periodicity cannot be substantiated, mainly because the 
random aspects of the structure of the polyphasic medium are very relevant, as might be expected in the 
case of polymer matrix in the active droplet. 
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(7.1) A summary of some results of averaged field theory 


In Figure 7.1 (a) a portion of a two-phase medium can be seen. In the droplet case, phase @ corre- 
sponds to the matrix and phase / to the interconnected pore space filled by a watery solution. 


/ Esup (5CP)} 


P+ #cG;dy—> SCP) 


Figure 7.1 (a) Figure 7.1 (b) 


The method of averaging over a volume is based on the possibility of defining three length scales: 
microscopic d , mesoscopic ¢, and macroscopic L , that verify d« 0, «L. After the averaging 


procedure the description of the fields is done in the macroscopic scale. 


In Figure 7.1 (b) one possible procedure to determine a characteristic length scale of the pore phase 
and the matrix phase is suggested (Scheidegger, 1963) (Suarez-Antola, 1984) (Sudrez-Antola, 1991). 
Four random selected points are shown, P, and P, belonging to the interior of pore space, P, and P, 
belonging to the interior of the matrix. 

For each point P belonging to the interior of a given phase, a spherical open neighbourhood 
E (C 50 ) of centre C and radius © is identified, such that E (C 30 ) is the open sphere of maximum 
radius included in the phase to which P belongs. 

Figure 7.1 (b) shows several examples for the pore space in the matrix of the droplet. 

A microscopic length scale d = sup {5 (2 )} can be defined for each phase. (Here microscopic is not 
used in the sense given to this word in physics). Our main interest is in the microscopic scale of pore 


space, d,. 


Now, to generate the mesoscopic length scale, for each point P we define a representative cubic vol- 


ume region B(P) as suggested in Figure 7.2 (a). This is called a representative volume element 
g eg g 


(RVE). In fact, for the purpose it could be a sphere or any other suitable region to apply the averaging 
procedure. 
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The point could belong to the interior of one of the two phases or could be located on the boundary 
between them. In the example shown, the point belongs to the matrix (phase @ ) but as already said, 
our main interest will be in the other phase. 


Figure 7.2 (a) Figure 7.2 (b) 


The length scale ¢ of this volume element is defined as follows. The RVE can be expressed as the 


union of two sets, each one corresponding to one of the phases: B Le ) =B Ae ) UB (P ) 


An indicator functi be defined for each phase: J,,(Q) Lf eB, (P) ( B) 
n indicator function can be defined for ea ase: = =a, 
c ction c ch p y 0 if OB, (P) M=a 
V(B_(P 
Then Hae La I ,(Q):dV = won = ¢,(P) is the volume fraction of phase jin the 


RVE B(P). 
Due to their definition, ¢, (ee ) + p, (P ) =] at each point of the polyphasic medium. 
Let yw be a scalar field in the active droplet (for example a mechanical stress that is a component of 


the stress tensor, or a hydrostatic pressure, or a concentration of a reactant). The intrinsic phase aver- 


ages of the field in the point P are defined this way, being y, (Q) = y(Q) ‘L, (Q) : 
1 
y,). (P)=——-~ y(Q)-1,(Q)-dv (u=a,f) 
( aA ) 7(B,(P)) a ( ) ak ) 
Figure 7.2 (b) suggests what could happen, when the length of the cubic RVE varies, when the intrin- 
sic phase average is (Vv, , (P) . If ¢ is too small and the point P does not belong to phase @ the 


average is zero. Then, as increases the average increases first and the oscillates until a length @, is 
attained, such that when ¢ is not too great compared to ¢,, the phase average (VV, : (P) remains 
almost constant. If ¢ increases even more, (WY. ), (P ) begins to vary relatively slowly, in a length 


scale L: when L>> ¢, >>d we define ¢, as the mesoscopic length scale and L as the macroscopic 


one. 
We refer all this to pore space. 
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The extrinsic phase average is defined relative to the whole RVE volume: 
1 
(v),(P)=(v,)(P)= VB Py dees” (2)tn(Q)-av =$,(P)(v,)(P) (u=a, Bp) 


If we average again first obtained average over the same RVE, we obtain the original average plus one 


f 2 
term of the same numerical order of : ) (Whitaker, 1999) that can be neglected. 


The original microscopic field can be decomposed: y= (y ) +W 
The extrinsic average of the gradient can be substituted by the gradient of the extrinsic average 
(V y )(P) 2V (y)(P) zO-V (y (P) (in the last relation appears the intrinsic average for a con- 


stant volume fraction ¢ of the considered phase). 


The Fick-like law J = —D-Vy can be averaged and approximated by (J ) = —D,y eV (y ) where 


Duy is an effective mesoscopic property resulting from the microscopic property D after applying a 


closure operation (Scheidegger, 1963) (Suarez-Antola, 1984) (Suarez-Antola, 1991) (Whitaker, 1999) 
(Davit et al, 2013) (Quintard, 2015). 
Even if the porous space behaves as isotropic at the microscopic scale, when the averaging process is 


applied, the resultant effective property De may be a tensor quantity if there is a degree of anisotropy 


related with the geometry of pore space. 

In the simplified models of the mass transport processes constructed in the following subsections, both 
in the droplet and the in the environment, we suppose that the effective properties are scalar fields at 
the mesoscopic level. 


(7.2) A model of the average elongation of an active droplet 


Let us suppose that a slight perturbation deforms the initially spherical droplet of radius 7, into a pro- 
late spheroid (a prolate ellipsoid of revolution, being z the axis of revolution) of major axis 


% 


Vl+eé 


a=h: (1 + é) and minor axes b= so that the volume V of the droplet does not change during 


the deformation: 


Vaden =a rab’ (7.1) 


Our purpose now is to describe the change in the semi-major axis of the spheroid applying the formu- 
lae (A.3.12) (A3.13) and (A3.14), that appear in Appendix 3, for the relation between the average rate 
of strain with volume and surfaces forces. 


We assume that the main average volume forces are the so-called filtration pressures p y acting on the 
polymer matrix (Scheidegger, 1963): 

By =V(9-P) (7.2) 
Here @¢ is the porosity of the droplet matrix and pis the intrinsic averaged pressure of the multicom- 


ponent solution in the pores. 

The local pressure in pore space, before averaging, can be approximated by a multivariable polynomi- 
al in terms of the concentrations of the components of the solution (including the solvent) beginning 
with linear terms (Green, 1969) (Suarez-Antola, 1997) (With, 2013). 
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When @¢ grows and tends to 1, the matrix disappears and all the drop interior is a multicomponent 


watery solution, possibly enclosed by a membrane. In this case the model is not applicable, and an ab 
initio viscous fluid dynamics approach to the resultant reaction-diffusion system is necessary. 
If 6=0 the pore space disappears, and the model is inapplicable again. 


So, we consider porosities ¢ not too close to one or too close to zero. 


The filtration pressure originates in a scalar potential field y =@-p (formula (A3.13) and the contri- 
bution of the volume forces due to filtration pressures can be calculated by formula (A3.14). 
We further assume that the main surface forces are due to surface tensions 7, , that could vary due to 
varying concentrations of diffusing reactants. 

1/1 | eee . .) 
If H= es a + ra is the local mean curvature (here R, and R, are the principal radii of curva- 

1 2 

ture) the second integral in (A3.12) can be written, despising, for the moment, Marangoni's shear 
stresses (Gallaire and Brun, 2017): 


Lon, Lfen,ern,) pare (058d (odene yg) att) 


0 


(7.3) 


92% 


0 
Then from (A3.12), approximating the average rate of strain deformation (2 é by the relative 
t 


eas sgh. 
rate of variation of the major axis of the spheroid —- ee a we obtain: 
a at 


ld 1 1 
eee a rer -@ en——- -e@ en -e@ en). -D ; 4 ‘ 7.4 
Parr 3-9 -V(B) fale Eon, (x é.ent+y-é, ’)| (¢-p+2-H 1) is} (7.4) 


Now we are going to establish relations between a certain reactant concentration and both, the aver- 


aged pressure p and surface tension v,. Through the chosen reactant, the distribution of both p and 


y, can be related with a measure of droplet chemical activity, and as we will see, to droplet stability. 


(7.3) A simplified reaction-diffusion model 


As we said before in section (4.1), there is a simplified approach introduced by Rashevsky to find an 
approximate solution to some diffusion problems, making some drastic approximations. Although in 
the case of a spheroid, some exact analytical work is feasible and has been done (Young, 1939 (a) and 
(b)) (Rashevsky, 1960)) we consider a droplet of roughly oblong shape with half-width r, and half- 


length r, (In our spheroid r ,=band r,=a). 


All local concentrations, transport parameters and rates of reaction in the pore space of the droplet are 
considered in the sense of the averaged field theory. 

As consequence, from now on, when we use the words “mean” or “average” in relation with these 
local concentrations of reactants in the droplet, we will be referring to the process that begins with 
already extrinsic average values at the mesoscopic scale and takes global averages at the macroscopic 
scale, sometimes on the whole droplet volume. 
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The previous averaging at the mesoscopic scale allowed us to homogenize an initial two spatially dis- 
joint phases in the initial mathematical model of the droplet. 

The permeabilities and diffusion coefficients are effective scalar properties as mentioned in the sub- 
section (7.1) above. 


Let the mean concentration of a certain reactant, roughly halfway between the centre and the periph- 
ery, bec. 


Besides, let the average peripheral concentration inside the deformed droplet at the ends be c, and at 


the sides be c,. The corresponding values just outside the droplet boundary are G and G,, respective- 
ly. 

We assume that across a distance 0 , measured perpendicularly to the droplet boundary, the concentra- 
tions e and or vary and reach a limit value c, corresponding to the bulk of the environment. 

We suppose that 6 is f the same order of numerical magnitude that the droplet dimensions (for exam- 
ple 0 =F, ). 

Assuming that there is no accumulatio or depletion of the reactant at the droplet boundary the follow- 


ing equalities can be posed: 
For the ends: 


LD pn \ 2:D, (4 
(2-4) =P-(¢ c')= 1G c.] (7.5 a) 
i, O 
For the sides: 
2-D, = ’ 2-D, ' 
L.(7-c)=P-(c,-c) = (c, -<,] (7.5 b) 
h O 


In equations (7.5) cross diffusion due to other components of the solution effects are neglected. 

The permeability of the boundary is represented by P (the boundary could a kind of membrane) and 
the averaged diffusion coefficient of the reactant inside the oblong droplet is D, . The diffusion occurs 
in the interconnected space of pores but, as already said, the system will be described using averaging 
field theory as a spatial superposition of the matrix with the space of pores. 

Due to this, D, is not simply a molecular diffusion coefficient, is an effective parameter that appears 
after fulfilling a closure process (Whitaker, 1999). 

The diffusion coefficient in the environment is D, . 

Equations (7.5) are a generalization of equation (4.1.3). This last equation was posed for a spherical 


droplet of radius €=7, when the diffusion in the environment is fast enough relative to the diffusion 


inside. If this happens also in the oblong droplet, then C, = G =c, and the equalities (7.5) at the ob- 


long droplet boundary reduce to: 


BE esp (ae) (7.6 a) 


si 
ETN =P (6) (7.6 b) 


N 


ty 
When the transport through the interface is fast enough in comparison with the diffusive transport in 


the environment and in the interior of the droplet c, = on and c, = e SO: 
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pag =a c,) (7.7 a) 
2D ps 2-D 
Brey zh) (7.7) 


The area of the two ends of the oblong droplet can be approximated by 2- (x . r;) and the area of the 


sides can be approximated by (2-4-r,)-(2-7). If the droplet geometric dimensions vary slowly 


enough relative to the reaction-diffusion processes, and if Q is the net rate of reaction producing or 


consuming the considered substance, the balance equation for the reactant can be approximated by: 


4 d_ (4 24D... ge 2D ys 
[Farin Se a( Sarin 0-2(w1f) Pee) (2-4-5)-(2-n) —P(e-e) 
This last mass balance equation can be recast as follows: 

$2 -9-3-7,| C9). E) (7.8) 
dt i ip 


If the reactant is the nutrient in the anisocoric generalization of Bertalannfy’s model of open system, 
c=c,, c,=X and Q= —(k, +k, ) -C, +k,-Cg Where Cc, is the assumed constant concentration of 
constituent, in this case the loose polymeric matrix in the droplet. 

If the reactant is the excreta in this model, c =c, , c, = Z and Q=k,-c, 

From now on, but only for the remainder of this section of the article, we will assume that QO is known 


and remains constant. 


In general, from boundary conditions (8.5) both, the interior concentrations at the ends (at the poles of 
the prolate spheroid) c, , and at the sides (the equator of the prolate spheroid) c, , can be expressed 
as functions of the mean concentration c inside and the concentration c, in the environment far 


enough from droplet’s boundary. Then equation (7.8) reduces to: 
d 2g (¢-&,) 
—¢ =Q-~—+ (7.9) 
dt A 


The characteristic time A is given by the formula, in terms of dimensionless numbers: 


a=( et _——.scS (7.10) 
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For a prolate spheroid of semi-major axis 7, = 1%: (1 + é) and semi-minor axis 1, = ——2— the charac- 
Vi+e 
teristic time is a function of the dimensionless deformation parameter € : A=A (e) 


Now, if we suppose that this characteristic time A(é) is at least an order of magnitude smaller that 


the time scale of the initial deformation of the droplet, the reactant will be always near its steady state: 
C(e)=c,+A(e)-O (7.11) 


In steady state the concentrations at the poles c, and the equator c, are given by: 


c(e)=c, +w,(e)-A(e)-O (7.12 a) 
c,(€)=c, +w,(e)-A(e)-O (7.12 b) 
In these formulae: 
21082) 
w,(€) = ——_—+—¥__ (7.13 a) 
: 21082), Pal6) 
D, D, 
282) 
w (ae (7.13 b) 
218?) Pale) 


Two cases are of interest for us. In the first case the resistance of the interface to the transport of the 
reactant is negligible in comparison with the resistance related with the diffusion in the droplet interior 
and in its environment. In the second case the resistance of the environment to reactant diffusion is 


negligible but the resistance to mass transport, both in the interface and in the interior of the droplet 
can be limit steps. 


When the mass transport resistance of the interface between the droplet and its environment is negligi- 
ble, P — +00. Then, we obtain: 


-l1 


es Naa RGM NAG ee 


w,(é) (24 (7.14 b) 
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2 | 
D, 
ee el 2 eee (7.14 c) 


-1 


a(e)-(2} et (7.15 a) 
6 5 re) , ro (¢) orlé) Mi (€) 
P D, P D, 
2 
Ww, é) =————_, (7.15 b) 
OS 
D, 
Zz 
WwW é)= (7.15 c) 
2) 9 Ponte) 
D 


(7.4) Droplet stability to small isovolumic shape perturbations: critical radii 


Applying the same simplified approach introduced to find an approximate solution to the reactant re- 
action-diffusion problem to formula (7.4), that gives the relative rate of elongation of the prolate sphe- 
roid, we obtain the approximate expression: 


ld dé 1 


Tne de Tg |P (PIG()-Ple(2)))-| («)) | (4 (a(2))-2-7,(a(e))) 


ie 0) 


(7.16) 


rs 
Here a= r-(1 + é) =1, (e) (pea ss (€) . The porosity fraction @ of the droplet matrix is 


Vl+eé 


constant by hypothesis. 
In (7.16) we suppose that both, the local average pressure in the pores of the droplet matrix and the 
local surface tension, are functions of the peripheral internal reactant concentration (in the poles are 


functions of c, (€) given by (7.12 a), and in the equator are functions of c, (€) , given by (7.12 b)). 


When the deformation parameter ¢=0, the spheroid reduces to a sphere of radius 
i, (0) =% (0) =r,. The concentrations in the poles and the equator are equal c, (0) =C, (0) =Cy. 
Here: 


c, =c, +w(0)-A(0)-O (7.17 a) 
With 6 ~ 7,, from (8.10) and (8.13) we obtain: 
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A(0)=[ "0 Jo(u et) to (7.17) 
12-P D,}) D, 


21-4?) 
D. 
w(0)= (GATS) 
7 eee P(e 
D, D, 
So, from (7.17 a) (7.17 b) and (7.17 c): 
Wt Fee 7 
Co =C,+>:| —+— |: 717d 
S e 6 [4 | Q ( ) 


The average pressure in pore space, and the surface tension take the values p (@ ) and 7, (és ) : 


Developing the right-hand member of the approximate formula (7.16) (that gives the relative rate of 
elongation of the prolate spheroid) in powers of the deformation parameter and retaining explicitly the 
linear approximation only: 


a eo) (7.18) 


The sign of the coefficient 2, determines the stability of the spherical droplet under small defor- 


mations. These deformations give prolate spheroids of: 
Semi-major axis: 
4 (é)=a=n-(I+e) (7.19 a) 


Semi-minor axes: 


b 
nle)= een (13) (7.19 b) 


The development concentrations at the poles and at the equator give: 
G(8)=¢5+G(0)-s+0(2) and e(2)=¢,+-04(0)-8-+0(2) respectively. 
é é 


The development of pressure gives, at the poles: 


a _ d, d 

B(«,(€))=Ples) +5 (cs) 4 (0) +0(2) (7.20 a) 
And at the equator: 

P(62(2)) = P(e.) +2 (¢,) <2, (0)-8+0(2) (7.20 b) 


The analogous development of the surface tension gives, at the poles: 


dy, (és) d 


r(a(2))=7.(¢s)+—S 5-4, (0)-e+0(¢) (7.21 a) 
And at the equator: 
dy.(cs) d 
y,(c:(€))= 7.(c,)+ LAE). 4 o,(0)-¢+0(6) (7.21 b) 
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The linear stability coefficient is, by definition: 


ap cee e ly (a(c,(e))-B(e,(e)))-|Z (c,(e)) (%(a())-2-7. (a (e))) 
' 2m dé : 74,(€) r, (€) 
é=0 
(7.22) 
From (7.16), (7.17) and (7.18) and (7.19) if follows: 
1 d l(d d d 3 
es ial ee ase cee NRE CH ema omy NR gee 
"Qn ct a 1) ry Ee ‘i )) E i) Gel )} 4-7, ) 
(7.23 a) 
If the surface tension remains constant, formula (8.20) reduces to this one: 
1 d d d 3 
A.=-— '|¢:|—P -| —ce,(0)—-——c, (0 a 7.23 b 
= fo (Srce)}(La0-Le0}+ 2-7, | (7236) 
In this case the effect of surface tension is always a stabilizing one, and the only destabilizing effect 
could stem from the term [4 P(cs (4 C, (0) _ #c, (0) if it could be negative and in absolute 


value greater than 


Ys 


sf 


Let us consider that a solution of the reactant in water fills the pore space. The local pressure can be 
approximated by the following formula, where now c is the microscopic solute concentration and 


c,, 1s the microscopic solvent concentration (Green, 1969) (With, 2013): 


P(c.c,, = [R-T-c-l,, -c-T,, ec, || R-T +c, - Ly °C, —Tyy Ce | (7.24 a) 
In this formula J, =J,,, so: 
& P(ese,)=RT=2-1,, 06-21, °6, (7.24 b) 
Cc 


Averaging (7.24 a) at the mesoscopic scale we obtain a local averaged pressure in terms of local in- 
trinsic averaged concentrations: 


Pel RE (eat Ao), Bele) (Ge), || RPG) BAe NC. tec), | 


From the molecular theory of liquids, we would expect that P(cs ) <0, due to osmotic effects 


Gs 
d(c), 


in the space of pores at the conditions of mild temperatures assumed to be prevailing in the model 
(Green, 1969) (Sudrez-Antola, 1997) (With, 2013). 

So, we assume this last inequality and return to the use c as an extrinsic phase average, as we did 
previously. 


d 
If the reactant acts as a surfactant, that is if ER Y, (ce ) <0 and besides if (8.23 a) is verified, then the 
C 


effect of the spatial variation in surface tension opposes the effect of the variation of pore pressure in 
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the periphery of the prolate spheroidal droplet and tends to stabilize the original droplet shape (a 


sphere). 
From formulae (7.12): 
#4 x(6)- [au (6)-A(e)| 0 (7.25 a) 
dee de" aa 
d d 
Ae (0)=[Zv(2)-a(e)) -O (7.25 b) 


The threshold of linear stability is attained when 4 ,=0 and instability of the original shape of the 
droplet is assured when 1,>0. 
When surface tension remains constant, from (7.21) and (7.25) the threshold condition 1,=0 is 


equivalent to: 


d d 3 
“(wy A) -—(w,-A)_ =| —— || —» —_ (7.26) 
dé é=0 dé é=0 4-7r,-Q b-|2 Bc ) 
dc 7 
If the surface tension varies between the spheroid poles and the equator and the increase in the local 


site ot , d 
reactant concentration in the drop’s periphery lowers the surface tension (— 7, (cs )). A ,=0 is at- 


dc 


tained when: 


Aer cae cars) s[eral-(Ee) 


% 


When the mass transport resistance of the interface is negligible, taking 0 ~ 7, we obtain: 


fw ‘A) = 2/_ % (7.28 a) 
de’ ' /’=0 Q9\D +2-D, 
d 1 r 
mar Peed ee ee 7.28 b 
ae” hea 9 sig | ey 
Then: 
d d 1 r 
anc A peck A ies 2 ON 7.28 
fal Mes Sela 3 aa | 280 


When the mass transport resistance of the environment is negligible, with 6 ~ 7, we obtain: 


%P 


d I D, 
— “A = peels : 8.29 
rus ) (=) n-P ( a) 
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ty P 
# (w,-A) _(_% \|_ 2 (7.29 b) 
dé =0 9-P hiP i, 
Then: 
ae oe 
= (wi), == (we a), =-|—2 || (7.29.6) 
dé 0 de a 9-P)| HP 


Considering formulae (7.25 a), (7.25 b), (7.28 a) and (7.28 b) the sense of variation of the concentra- 
tions in the poles and in the equator is the same in both cases, when P-—>+0o and when 


D, — +00, and depends of the sign of the net rate of reaction Q producing or consuming the con- 


sidered substance. 


When Q>0O we have 4 (0) = (2 w (é) . A(e)) -Q<0O (the concentration decreases rela- 
é é e=0 
d d Ne aod 
tive to c, in the poles) while ae ee (0) = We w,(é) “A () - Q>O (the concentration increases 
é é é=0 


relative to c, in the equator). 


d ‘ ; ; . 
When QO <Owe have ie (0) >O (the concentration now increases relative to c, in the poles) 
é 


ae: ' 
while de Cy (0) < 0 (the concentration decreases relative to c, in the equator). 
é 


When the resistance of the interface to mass transport is negligible (P — +00), and the net rate of 
reaction Q is negative (the reactant is consumed inside the droplet), from equations (7.28 c) and 
(7.27) we find: 


. 9-(D, +2-D,) 
(7.30) 


il eae’, ae raid ee vane 


£ r (es) 


d d 
If ae a ic )is always zero or can be neglected relative to aE P(cs ) , from (8.30) we obtain the fol- 
Cc C 


lowing formula for a threshold radius: 


1 
-|2] (D,+2-D,)-y, (731) 
y) 


Once the droplet has grown to this radius and surpass it, the spherical shape gets instable, and the drop 
begins to deform progressively adopting the shape of a prolate spheroid that elongates without chang- 
ing its volume. 
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A comparison between formulae (7.30) and (7.31) shows that a variable surface tension that decreases 
when the concentration increases, increases (all other conditions being equal) the threshold radius. 


When resistance of the environment to mass transport is negligible the equation that allows to deter- 
mine a threshold radius, is this one: 


a 27-D, Vs (cs) 
—3 -|- , } eee (7.32) 
ea FI pelt ene) 


} lent) 


P(cs ) can be neglected relative to the contribution 


d= 
Up to now, we have assumed that @- (|¢ Dp (e ) 
Cc 


If the contribution of filtration pressure ¢ { 


dc 


ld 
of a variable surface tension —-—y, (a ) , the linear stability coefficient (7.24 a) is now given by: 


% 


an-elf-+ fred} (Za-La()}-F-r(e)| (733) 


0 


Nevertheless, a variable shear tension is associated with variable peripheral shear stresses (Marangoni 
tangential stresses) whose effect on the relative rate of elongation of the prolate spheroid is not consid- 
ered in (7.33). 


To do this by analytical methods, let us consider the solution c.(0 ) of the steady reaction-diffusion 


Q 


equation yet =0 in the surface of a prolate spheroid, when the resistance of the interface to 
i 
mass transport is negligible (P—-+0o) and the concentration c, in the environment verifies 


V°c, =0. Here, as usual, V (2, =VevV ee is the Laplace’s operator. 


The angle @ is that formed between the z axis direction (z is the axis of revolution of the spheroid) and 
the straight-line segment that connects the centre of the spheroid with a point of its surface. 


The boundary conditions for the partial differential equations are the continuity of the concentrations 
and the continuity of the normal components of diffusion flows, when we cross from the environment 


to the interior of the spheroid, and the presence of a limit concentration for c, when the distance to the 


spheroid increases without bounds. 


With these boundary conditions we can adapt the solution originally derived by Gail Young (Young, 
1939) as follows, in the first order in the dimensionless parameter ¢ : 
Qn 


-cos’6 7.34 
(2-D, +3-D,) — oe 


c,(0) ¥ Cs -E- 


If the surface tension is given by the formula /, (c.) =/, (cs ) + B- (e, —€ ) , from (7.34) we obtain: 
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2 
0-1, 


————$——-: °@ 7.35 
6p30) ee 


(C=, (€s)-E-B: 


When £ is negative, if the net rate of reaction Q is positive, the surface tension has its maximum at 
the equator and its minimum at the poles of the prolate spheroid. If Q is negative the surface tension 


has its maximum at the poles and its minimum at the equator. 


From (7.4) we can separate the contribution of the surface tension to the relative elongation rate for- 
mula. Considering the tangential stresses acting on the spheroid surface when the surface tension var- 
ies in this surface, we can adapt again a solution constructed by ae (Young, 1939) to obtain: 


d 1 Ny 
een Qe ft .G4jo pe — | 1.36 
a), uf 7.(¢s)- r, ce OD = D 5) epOKe) 2) 


At the first order in €, the contribution Wage of both, the variation in surface tension due to the 
t 


combined effect of the variation in the surface concentration of the reactant and the transmission of 
(ae Om 
n (2-D,+3-D,) 


When / is negative, if Q is negative the corresponding contribution to the relative elongation rate is 


Marangoni’s stresses, is — 


negative. So, considering only surface tension related effects, small shape perturbations cannot desta- 
bilize the originally spherical shape of the droplet. 
But if we have a positive net rate of reaction Q, and / is negative, then we have a positive contribu- 


tion to the relative elongation rate for a small shape perturbation. 
In this case there is a threshold radius to linear instability due only to surface tension effects: 


¥,(cs) (2:D,+3-D,) 
(-4) Q 
The solution (7.34) for the concentration field in the boundary of the droplet, and the equation (7.24 b) 
for the derivative of the pressure as a function of this concentration, allows us to calculate a linear 


(7.37) 


"th = 


approximation in terms of € to the variation of pressure in the inner side of the boundary. 

Then, the approximate (for ¢ <1) formula Ca a P(c;): “ce, (9)-e+o/(e) can be intro- 
duced in formula (A.3.14) for the effect of volume forces that depend of a potential field y, in this 
case being y =@- [> (cs)+¢- “> (5) ' c (0)-e+0(8)] ; 


The integration on the closed boundary of the assumed constant term @- P(cs) is zero, SO we obtain 


for the relative elongation rate due to volume forces: 


al x Satay ll [e227 F: (re 049-4 07)}(6-£ wl.) Selo) ash or o(e) 
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After calculating the surface integral, the result can be added to formula (8.36) and a more rigorous 


formula for the relative elongation is é= @ é| + a: € | can be derived. 
dt dt ), \dt ), 


We have seen in this section that, once a threshold is reached, the spherical shape of the drop gets in- 
stable and a progressive isovolumic (isochoric) elongation occurs. We obtained closed form analytical 
formulae for the critical radius with which the stability is lost. 

When the dimensionless parameter € , that gives a measure of the isochoric change of shape in the 
droplet, is not any more little relative to 1, the linearized mathematical model discussed in this section 
should be expanded to account for the flow of the solution in the pore space, how the pore geometry 
changes during deformation of the droplet and how its viscous matrix deforms and flows. 


A straightforward development of the member to the right of the equality sign in formula (7.16) or in 
the less approximated formula obtained adding (7.36) and (7.38), in powers of the parameter €, to 


nell sets 
give oH E=A,-et+d,-€' +A,-E>+... would not be useful in this case. 
t 


The reason is that the fluid dynamics aspects have not been explicitly covered in this section. 

The full coupling of fluid dynamics and chemical reactions (including electrochemical aspects) 
needs to be included to completely describe the active droplet behavior, including droplet division, 
droplet movement and droplet fusion. 


Furthermore, cooperative effects are related to the self-assembly of a droplet germ from a solution and 
the eventual development of sequences of spontaneous pattern formation inside the active droplet dur- 
ing its growth. 

At critical values of the parameters, when nonequilibrium phase transitions occur, and in chemical 
networks in which the number of certain molecules is small enough, random effects must be included 
in the mathematical models, following the guidelines summarized in section (9.2) of Appendix 1: 
“Random dynamic systems”. 

These problems of biological interest will be considered in the next section, dedicated to discussion 
and final conclusions. 


(8) Discussion and Final Conclusions 


(8.1) The model of open system with a single autocatalytic hemi-reaction proposed by van der 
Vaart. 


Let us begin with a historical comment about autocatalytic reactions in growth equations that are relat- 
ed with the results obtained in section 3 of this paper. 

We have seen there that the analysis of the stability of the steady state point of the open system model 
with a single autocatalytic hemi-reaction, proposed by van der Vaart, leads in all cases to a stable fo- 
cus. 

So, the product of the reaction cannot present the sigmoid behavior sought by supporters of the auto- 
catalytic limiting step as a determinant of growth kinetics, as they did it mainly since the early twen- 
ties to the thirties of the last siécle. 

If open systems had been used from the beginning, the autocatalytic model would not have acquired 
the importance it had in the history of the global mathematical description of growth and in population 
dynamics. 
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However, autocatalytic kinetics are essential to understand the details of the dynamics of the delicate 
biochemical machinery at the organism level, and the complex interaction between genotypes, pheno- 
types, and environment at the population level. 


(8.2) Generalizations of von Bertalanffy open system and a study of its anisocoric versions. 


Now, we reconsider first some of the results obtained in section 4, subsection (4.1), (4.2), (4.3) and 
(4.4). 

Then, we reconsider the results obtained in sections 5 and 6. There, the kinetics of open anisocoric 
complex reaction-diffusion systems with a limit step to mass transport located at their boundaries is 
formulated, and slow manifold theory is used to derive von Bertalanffy’s growth equation for these 
systems. 

Some properties of stochastic dynamic systems are briefly considered. 


(8.2.1) Generalization of von Bertalanffy's model to include variable volumes and bulk diffusion 
with a nonconstant effective permeability 


In subsection (4.1) an effective permeability P, was introduced to combine the effect of the bulk dif- 


fusion of a chemical species inside the system with its permeability due to the boundary barrier: 


Py = 


Here D is a diffusion coefficient of the substance, ¢ is a characteristic length (perhaps the radius of a 
sphere) and P is the permeability of the boundary to the considered chemical species. 
This effective permeability, if it turns out to be applicable, allows to avoid the complications associat- 
ed with the use of partial differential equations to describe reaction-diffusion systems. 


When the dimensionless number 


is small relative to 1, the limit step to mass transport is locat- 


ed in the boundary and PR, =P. When is at least and order of numerical magnitude greater 


2°D 


than 1, bulk diffusion dominates mass transport and RP, ® 


When the system consists of a population of dividing spherules (a biophase), as shown in () it is possi- 


ble to choose a representative value of the characteristic length (=, and work with an effective 


oe 
[1 
2-D 


The effective permeability appears in the critical value X. of the external concentration of nutrient for 


permeability P, = that do not changes during population growth. 


the open anisocoric system without scale effect, so is a main determinant of the bifurcation that sepa- 
rates a self-consuming from a growing and dividing system. 


What happens when the system is not able to divide, and the volume can change at least one order of 
numerical magnitude during growth? 
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The simplest geometry is a sphere. So, let us consider an open spherical system of radius 17, (t) that 


has the same structure of kinetic equations as the anisocoric generalization of Bertalanffy’s model with 


scale effects. 
The diffusion in the environment is fast enough and external source of nutrient is big enough to have a 


constant external concentration X . The now spherical volume is proportional to the number of moles 


1 2 
of constituent V =——-n, so if C, is an average value of nutrient concentration the equations of evo- 


Ce 
lution are: 
d 
ret S +P; -(X -G)+(-(k, +&)-@ +k, -¢5)-V (8.1) 
d 1 
—V =—(k,-¢,-k,-cg):V (8.2) 
dt Ce 
Now: 
4 
VS et (8.3 a) 
S=4-1-r (8.3 b) 
P 
P; (1%) =————~ (8.3 c) 
14% ud 
2-D 
If the scale of time of growth is small enough relative to the scale of time of variation of the nutrient 
concentration: 
d_ § = = 
qivy fe (X-G)+(-(h +e) +h -¢o) (8.4) 


In the time scale of volume variation c, can be considered in quasi-steady state. From (8.4): 


3 
= Po( 1) X +ky ty Pa (to): X + hy Co 
C, =-9 (8.5) 


cia S 3 
ky +h + Po(r) k, +k, +—-Py(%) 


% 


Substituting (8.5) in (8.2) and after some reordering, we obtain the following equation for a von Ber- 
talanffy’s the growing droplet: 
3 
—-P,(%)-(X -X,)—-k, +k, Cg 
d SA lcr, 
Ge) a (8.6) 
eS k, +k, +—-P (1%) 
0 


In (8.6), the effective permeability P, (%) is the function of droplet radius given by formula (8.3 c) 


k 
and X, =—: C@ is the same critical concentration that we found in (4.2.2) “Scale effects in the aniso- 
1 
coric generalization of von Bertalanffy’s model of open system”. 
Like Equation (4.4.9) for the volume of the system, Equation (8.6) for the radius of the droplet is an 
external approximation in the sense of singular perturbation theory (see subsection (9.3.3) “Slow man- 


ifolds” in Appendix 1. 
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The steady state radius (fixed point, rest point or equilibrium point) of (8.5) are these two: 7% =O for 


2) 6-k,-D 
+ 


D 
any value of X and %y =,/| = (X =X. ‘J —— with physical sense if and only if 
P KeaKy P 


X2X,. 
If X < X, the rest point % =O is globally asymptotically stable and the system self-consumes. 


If X > X, the rest point % =O is unstable and the rest point is stable: 


This rest point attracts all trajectories with physical meaning whose initial point 7, (0) verifies 
ty (0) # (). If the initial radius is small enough, the system can grow tending asymptotically towards a 


maximum radius 7, y, - 
The external concentration of nutrient is a control parameter also in this model that we call “von Ber- 


talanffy’s growing droplet”, and there is a bifurcation when X = X,. 


(8.2.2) Self-consumption and growth in generalized anisocoric models of open physical-chemical 
systems with and without scale effects. 


In subsection (4.2) the original model of open system due to von Bertalanffy, including a nutrient, a 
constituent and one excreta, was generalized to include a volume dependent on the composition of the 
system. 


; n 
The volume of the system was taken proportional to the number of moles of constituent: V=—. 
Ce 


This implies a constant concentration of constituent c, and the independence of the system volume of 
the number of moles of nutrient and excreta. 


2 3 
When the scale effect was included, so S = u-V° (0<6 <1, for example 6 =— or 6 =—) and the 
H p 3 A 


external concentration X of nutrient was taken as control parameter, the dynamical system describ- 


ing the time evolution presents a bifurcation for a critical value X, =—*-c, of the control parameter. 
1 


If X < X, the system asymptotically disappears, but when X > X, it tends to a limit steady state 


1 
whose volume is given by V,. aut BEE fixe Soe) ' Thus, for ne we have 
We ne ge || oie 3 
3 
"Ce | ky ok, ; 


If the kinetic constant k, of the catabolic reaction tends to zero, the steady state volume V,, grows 


without bounds. So, according to this mathematical model, the catabolic reaction is a key component 
in growth control when the system, considered as a highly simplified model of a single organism, in- 
evitably presents a scale effect. 

Numerical simulation of the dynamics suggest that sigmoidal growth of system volume is produced 
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when the initial volume is small relative to the steady state volume. 

Besides, numerical simulation ° suggests that when the initial number of moles of the system’s com- 
ponents are small in comparison with their steady state values, their evolution follows two-time scales. 
One fast, that appears first, and the other slow that dominates the asymptotic approximation to the 
steady state. 

During this second stage the concentrations of nutrient and excreta are nearly constant and equal to 
their steady state values, so the mass of the system can increase as a function of time, following a sig- 
moid curve also. 

Qualitative methods for phase plane were applied to study the geometry of the orbits 
when X < X,and when X > X,. The corresponding pattern of orbits, including the separatrices, is 
shown in Figures 4.7 (a) and 4.7 (b). 

An application of the theory of topological indices in the plane (Roxin,1972) (Attle-Jackson, 1989) 
(Jordan and Smith, 2007) provides qualitative information on which bifurcations are possible by vary- 
ing a control parameter, in this case the external nutrient concentration X . 

It is possible to define a partial Poincare’s indices for a rest point of a plane smooth vector field 


f (x ie), working with the orbits located inside sectors bounded by separatrices that end or begin in 
the rest point. 

A conservation of the sum of partial indices of a set of rest points in a suitable subset of the phase 
plane, bounded by separatrices, can be employed to predict possible bifurcations and stability proper- 


ties of fixed points. This can be done in the case of the plane dynamical system obtained from the an- 
isocoric generalization of Bertalanffy’s model, with scale effect. 


When the scale effect is absent, S = -V and the dynamics predicted by the resultant mathematical 


model is like that predicted by the mathematical model of biophase due to Perret and Levey, but cer- 
tain interesting similitudes and differences between these models appear. 

Both, the biophase model of Perret and Levey and the biophase model developed in the present paper, 
are anisocoric models with no scale effect, which can be interpreted as a population of spherules that 
divide asynchronously. 


The model of open and anisocoric system due to Perret and Levey is based on a single chain of ana- 


bolic first order reactions with no branches A,z2 A, &...<2 A, , G2A, . The volume of the sys- 


: ite : 1 
tem is proportional to the number of moles of the last chemical species in the chain V = —-n “ss 
Ce 


Only A, can cross the boundary of the system: the other components are constituents. The external 
concentration X of A, is assumed to be constant and is a control parameter. 

As summarized in subsection (9.1.3) of Appendix 1, the kinetic matrix of the linear and non- 
homogeneous system of differential equations Perret and Levey K (Xx ) is a matrix of a class known 


as Metzler matrices. These matrices have a dominant real eigen value A of x ) such that the real parts 


of the other eigenvalues are strictly smaller than / o( X ) . When the external concentration X of A, 


is greater than a critical value, the asymptotic growth of the system volume is exponential with coeffi- 


cient A(X). 


® Done in 1976 by Eng. Carlos Zamalvide Garderes and a couple of years later by Eng. Fernando Rodri- 
guez Bogorja. 
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Also, in the anisocoric generalization of Bertalanffy’s model of open system, without scale effects, 
studied in the present article, the external concentration of nutrient X can be taken as a control pa- 
rameter and when this concentration is greater than a threshold value the asymptotic growth of the 
system volume is exponential. 


However, besides the number of reactants, the main differences between the two models can be found 
in the catabolic reaction, that in Perret and Levey’s model do not appear. 

Let us then leave the Perret and Levey biophase model and focus now on the biophase model proposed 
in this work. 


k, 
MPs 


The corresponding dynamic system presents a bifurcation when X = X, =| 1+ -X, with 


* 


X, =—-C, being the critical external concentration of the anisocoric model with scale effects, and 
1 


k being the kinetic constant of the catabolic reaction. 
Thus, all other conditions being equal, and for k, #0, the critical concentration X, that must be ex- 


ceeded so that a system (or a single microorganism), which does not divide, can grow to a maximum 


size, is always lower than the critical concentration X, that must be overcome so that a biophase (or a 


population of microorganisms) can enter a mode of exponential growth. 


In the absence of catabolic reaction, that is, when k, =O it turns out X = X,. 


But whenk, #0 , X, depends on the size of the spherules indirectly through the parameter yu (the 


relation between yz and the volume of the spherules was found in subsection (4.3)), and directly of the 


a ee 
eas 
2-D, 


permeability of the membrane or interface and of bulk diffusion in the interior of the spherule with a 


effective permeability P, = of the barrier to the nutrient (that includes the effects of the 


representative characteristic length ¢, to make P, independent of time). 


The critical nutrient concentration X, for a biophase with an irreversible catabolic reaction decreases 


as spherule size and the effective permeability increase. 


If X <X, the system asymptotically disappears, but when X > X.. it tends to an asymptotic un- 


bounded growth regime. In this regime the volume increases exponentially as a function of time: 
V(t)eV, -exp| 2,(X)-1] 


We found that the growth coefficient when X > X. is given by the formula: 


2 
aoe He Fach x x,) (k, +k, +k, +uU- Ps) 
2 Ce 2 


If the kinetic constants are related with the absolute temperature T using suitable formulae obtained 
from the absolute reaction rate theory, a relation between the average time between generations 
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A(X) 
1977) that is found in experiments that relate the measured average time between generations and the 
measured absolute temperature (Rose, 1976). 


and Tis obtained that is qualitatively (Figure 8.1) the same relation (Sudrez-Antola, 
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Figue 8.1 Qualitative draft of the variation of the average time between generations with absolute 
temperature for the model of biophase with catabilic reaction. 


(8.2.3) A probabilistic foundation for the linear relationship between surface and volume in the 
biophase. 


In the subsection (4.3) the quotient between the total area and the total volume of N spherules forming 
a biophase (at a certain instant of time) was introduced as a random variable: 


we 1 <<. 
gD DS 
cant ae N 194 N 

Doe ey 

k=l N k=l 


Here appear the sums of the random areas S$, and of the random volumes ¥, of the spherules. 
Both s, and ¥, were considered as equally distributed random variables with finite variance, being 
their mean values 5 = E|s| and v=E [7] respectively. 
Then the strong law of large numbers was applied to conclude that, when N +o, {converges al- 
most certain toa number “= = 
Fr 
So, for a numerous enough population we can expect that $ = w-V 
As the scale effect occurs during the growth of each individually taken spherule, it was assumed that at 


each moment of the life cycle of that spherule, between its area S and its volume ¥ a non-linear rela- 
2 


tionship of the following type is established: § ~ £4, -¥? 


Considering the inequalities 5; <S, <5,, ¥, <0, Sv, (K=1,...,N) and the definition of “ , both 
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- ' sats oes Ce dsp yer ee 5; . 
i and its almost certain limit 42 were bounded in this way: oe Sus =? ; a sus ee 
f i f 
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Then v; = 2-v, and the formula § © 44-0? were applied to connect the area with volume at the same 


instant of time for the same spherule. From the inequalities, lower and upper bounds were obtained: 
= cH a 
V;3 Ho V;3 


This last constraint shows that the parameter “/ increases when the initial volume v, of the spherules 


decreases. 
As consequence the critical environmental concentration of nutrient for the biophase 


k 
X= [ + A - X, decreases when the initial volume of the spherules decreases. 
Mile 


Besides, the formula for the growth coefficient A(X ) of the biophase, shows that when v, decreases 


this growth coefficient increases. 


(8.2.4) External approximation and slow manifold theory applied to the anisocoric generaliza- 
tions of Bertalanffy’s model of open physical-chemical system. 


In subsection (4.4) “A singular perturbation approach to the dynamics: the volume and mass growth 


equations” we obtained Bertalanffy's growth equation “ =qa-M°—x-M as an external approxima- 
t 


tion (in the sense explained in subsection (9.3.3) “Slow manifolds” of Appendix 1) to the original dy- 
namics of the anisocoric generalization of Bertalanffy’s open system model. 

But as we already said, there is an important difference because in the obtained growth equation and 
the original Bertalanffy’s growth equation: in this last case the anabolic constant @ and the catabolic 
constant «K are empirical parameters, while in the equation obtained here they are given by well- 
defined analytical formulas in terms of the more basic parameters of the model of the physical- 


= ip Hef (X, X,) 


* 


chemical open system: a K=k, 


Here f (xX joe .) originates in a more general formula J, = f [x ; 2 for the flux density of nu- 


trient, so that it can present the kind of saturation effects that are often observed in microbial physiolo- 
gy. 


(8.2.5) External approximation and slow manifold theory in the general framework of chemical 
kinetics of open systems. 


In section 5 “Formal kinetic models for open homogeneous reaction-diffusion systems with com- 
position-dependent volume” we introduced a generalized framed to study the general formal kinetics 
model of anisocoric open systems with a limit step to mass transport located at their boundaries (or if 
bulk effects due to diffusion inside the system can be included by means of Rashevsky’s approximate 
method). 
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In section 6 “Heuristic derivation of an equation analogous to the von Bertalanffy growth equation” 
we generalized the approach by singular perturbations and slow manifolds to more complex open reac- 
tion-diffusion systems whose volumes depend of their composition, and obtained again von Ber- 


d 
talanffy growth equation ae =a-M° —x-M as an external approximation to their complex dy- 
t 


namics, with closed form analytical formulae for the anabolic and the catabolic terms: 


nen 


C 5 


, with 7 = >, ;-J; given as a weighted combination of all the nutrients flow densities 


icy 
> M.-8; 


iel', 


through the system boundary, and «= given as a weighted combination of all the effec- 


fe) 
tive densities of catabolic reaction rates of each excreta divided by system mass density. 
All these terms are calculated on the slow manifold corresponding to the dynamics (assumed to exist). 


Now, how such expressions like > M ;-J, and > M ;-&, could maintain a relative constant value 
icy icl', 

during growth, as we find by adjusting the parameters of the growth equation to experimental data? 

A possible answer to this question, in a physicochemical system with numerous components and nu- 

merous chemical reactions between them, could be as follows. Let the existence of a slow manifold in 

the state space be admitted. The trajectories of the system, possibly after a fast transient (a variation in 

the inner time scale), get a neighbourhood of this slow manifold and remain close to it. 

The slow manifold has certain slowly varying concentrations (varying in the outer time scale), in 

whose proximities the concentrations of the components of the system evolve, perhaps with small 

fluctuations. 

Then, it could be expected that the variations in concentrations associated with the chemical reaction 

network, including oscillations in the form of limit cycles and even chaotic trajectories of some com- 

ponents, when considered from the global perspective of the density and volume of the system, com- 

pensate each other so that they do not influence the dynamics of the volume or the total mass of the 

system. This heuristic approach seems to warrant a detailed investigation and connects with the next 

subject to be discussed. 


(8.2.6) Stochastic effects in reaction-diffusion systems 


At critical values of the parameters when nonequilibrium phase transitions occur, and in chemical 
networks in which the number of certain molecules is small enough, random effects must be consid- 
ered in the mathematical models of the corresponding systems. 

Near phase transitions the random nature of the microscopic world (now in the sense of physics, not in 
the sense of averaged field theory) must be included in the description of the dynamics due to the long 
range of the correlations between the interacting particles. 

In the chemical networks of cells or prebiotic systems, the numbers of certain key reactants vary be- 
tween numbers less than ten and a few thousand. Due to the probabilistic nature of the molecular reac- 
tion steps, the local numbers change randomly and with significant relative increments. 

In this work we suppose that the evolution of the reactant concentrations can be decomposed in the 
sum of a deterministic term that stems from usual formal kinetics, and a random noise term with a 
well-behaved probability density (that is, with zero mean value and well-defined finite variances) fol- 
lowing the guidelines summarized in subsection (9.2) of Appendix 1, “Random dynamic systems”’. 
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This is known as the linear noise approximation. When the linear noise approximation can’t be ap- 
plied, there are methods of data analysis that in principle allow us to distinguish between chaos, almost 
periodicity and random noise. 

Although random fluctuations can modify the bifurcation scenario, under conditions of applicability of 
the linear noise approximation, it seems that a deterministic analysis as was done previously can ex- 
plain the stability behavior at the level of mathematical modelling chosen here. 

This was assumed in the body of this research, but in the discussion in section (8.4) below, relating 
prebiotic systems and the problem of the origin of life, the consequences of including random phe- 
nomena in mathematical models will be reconsidered, and their importance from both the selection 
and evolution viewpoints will be stressed. 

This will be done discussing prebiotic systems that spontaneously originate spatial patterns, are able to 
modify the balance of mechanical forces, grow and then divide in into in general two appropriately 
structured offspring to in turn grow, last, and divide. A preliminary discussion is given in next subsec- 
tion (8.3) “Active, growing, and dividing droplets”. 

The processes that counteract and overcome, for as long as possible, the disorganizing effects of the 
inevitable random fluctuations from without and from within, must be controlled by some mecha- 
nisms. 

It is not enough for the prebiotic system to have a membrane bound continent ’ and a metabolic net- 
work, it seems necessary to have a carrier of information capable of replicating itself and constraining 
physical-chemical events, again playing its role in both systems resulting from the division: in short, in 
words of Howard Pattee it is necessary to have “the symbol in the matter”. 


(8.3) Active, growing, and dividing droplets 


In section 7 after a summary of the volume averaging method, we applied it to construct a mathemati- 
cal model of growing active droplets with porous matrices. The stability threshold of droplet’s original 
form was considered under varying filtration pressures and surface tensions. 

Approximate analytical formulae were derived for the droplet’s critical radius as a function of a reac- 
tion rate, surface tensions at the drop interface and transport parameters of a reactant in the drop and 
its environment. 

When the filtration pressure effect is dominant and the surface tension is uniform, the critical radius 
occurs when the rate of reaction is negative (the reactant is consumed inside the droplet) and is given 


by the following formula: 


9 3 (D, +2-D,)-7, 
oth = 4 aa) |= | 


If surface tension variation along the boundary of the droplet is dominant, if Marangoni tangential 
stresses are included in the analysis, and if the reactant acts as a surfactant (in the formula 


Y, (c) =/, (a ) +B: (c 5 cy) the parameter / is negative), then the critical radius occurs also when 


the reactant is produced (Q > ()) inside the droplet, and is given by the following formula: 


7 Or another structure with similar functions, if the “membrane first” hypothesis is discarded, as think the 
proponents of a “coacervate world first” as a path to the emergence of life. 
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In this model net the rate of reaction related with the considered reactant is assumed known. 
However, it is possible to construct a similar active growth model for the anisocoric Bertalanffy’s 
model, with a reaction rate for each diffusible component, nutrient, and excreta. 


Let us review now some theoretical and experimental results, related both to the stability of active 
droplets, following a historical course, although not exhaustive. The focus will be on open reaction- 
diffusion systems related to the topics covered in this work. Then we can study their connection with 
the model constructed in section 7 of this paper, “Active droplet growth and stability”. 


We begin with Young’s analytical solution corresponding to a mathematical model of an open reac- 
tion-diffusion system already mentioned in section 7.3 (Young, 1939). 

For an elongating prolate spheroid he obtained the concentration fields for any value of the semi-major 
a and the semi-minor b axes. 

From the knowledge of these fields and assuming a_ variable surface tension 


Y, (c, (0)) =7, (c,)+6 Y, *COS *@, he derived a formula for the relative rate of elongation of the 


semi-major axis a due to surface tension related effects (Young, 1939): 


ld 7, (cs) (<2) oy.) 1 
—s =— . . = s 7 — 8.7 
(+ = a) a, Gps ee + (Qty —a5,) male (8.7) 


As in section 7 above, the angle @ is that formed between the z axis direction (z is the axis of revolu- 


tion of the spheroid) and the straight-line segment that connects the centre of the spheroid with a point 
of its surface. 


When the ratio a, increases, the coefficient a@,, that measures the effect of the transmission to the 


interior of the body of Marangoni's tangential stresses decreases, while the coefficient @,_ that 


measures the effect on this relative elongation rate of the spatial variation in surface tension increases. 
When azb , ay >@s, and for a~2-b both effects are almost equal, and after that @), <@;,. 


So, the contribution of the second term in (8.7) to the relative elongation rate changes its sign. 
As consequence, a net rate of reaction Q the tends to stabilize the system for nearly spherical shapes 


tends to enhance the elongation in the case of prolate spheroids with a >2-b 


When the ratio af increases the positive coefficient a, decreases, from 1.6 when a xb to 1.5 when 


a~2-b, approaching 1.18 when a>>b. 


Rashevsky, Young, Landahl and other researchers (Rashevsky, 1960) developed a (perhaps premature) 
open system model of biological cell division based on a volume force and a variable surface tension. 
The volume force is due to the push that the concentration gradients (due to the diffusion of the reac- 
tants) exert on the particles of the continuous medium that constitutes the interior of the cell. 

There is no polymeric matrix as such: it is a cell model considered as a membrane-bound bag of reac- 
tants that constitute a viscous liquid solution in which reaction-diffusion processes occur. 

For a constant surface tension, Landahl (Rashevsky, 1960) rewrote the formula for the elongation of 
this kind of cell model (developed by Rashevsky and Young) as a third-degree polynomial: 
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e=A-e-Bee*—C-e" (8.8) 
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The coefficient A changes from negative to positive when the radius of the initially spherical cell 
attains a critical value, but in general the other coefficients are positive. 
So, equation (8.8) suggests that the growth of € is bounded. 


However, this result does not mean that the division is not possible. 
If a constriction develops progressively in the middle of the deformed membrane-bound bag, this can 


go on being compatible with a constant asymptotic value ¢,, of the dimensionless shape parameter, 


because the constriction could generate a dumbbell-like shape with two smaller spheres (relative to the 
original unperturbed one) linked by an ever-thinner bridge. 

As suggested in Figure 8.2, the points in the middle of bridge surface have principal radius of curva- 
ture of opposing signs, while the points in the surfaces of the spheres have principal radii of the same 
sign. 


Figure 8.2 (modified from Gallaire and Brun, 2017) 


However, one of the principal radii of curvature in the points of middle of the bridge (red arrow corre- 


sponding to R, in Figure 8.2) is much smaller than the other, is positive and smaller than the positive 


radii of curvature (blue arrow pointing to R,.) in the points located in the spherical surfaces of the 


dumbbell. 

If the surface tension is the same at each point of the boundary, assuming a quasi-static relation (La- 
place’s law) between the pressure difference across a point in the surface (inner pressure minus envi- 
ronmental pressure), the mean curvature there and surface tension, the inner pressure in the bridge will 
be greater that the inner pressure in the spheres. 

The resultant gradient of pressure (horizontal red to blue arrow on the axis of Figure 8.2) will tend to 
enhance the separation. 

This is quite analogous to what happens in the Rayleigh-Plateau instability. This instability appears in 
an initially cylindrical thread of liquid whose shape is suitably perturbed and eventually divides in 


droplets (Gallaire and Brun, 2017). The cylinder with radius R, , outlined with dashed lines in Figure 


8.1, corresponds to the above-mentioned cylindrical thread of liquid. 


Deformations and instabilities in liquid droplets have been experimentally investigated since the nine- 
teenth century (Hanczyc, 2014). 
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There has been both experimental and theoretical research in instabilities of viscous and chemically 
inactive liquid droplets immersed in a liquid environment, caused by variations in surface tension (the 
so-called Marangoni effects). 

Suitable variations in surface tension produced by the localized addition of a surfactant in droplet’s 
surface, induces a viscous flow inside the droplet. 

The flow redistributes the surfactant and modifies the gradients of surface tension until the original 
shape of the droplet gets unstable and eventually divide (Greenspan, 1977 and 1978). 


Figure 8.3 shows the flow that sustains a deformation process in a viscous liquid droplet due to a 
symmetric gradient of surface tension. 


Figure 8.3 A draft of the viscous flow inside an inactive droplet due to Marangoni effects. 


A bulk flow of liquid is produced in the droplet interior that moves toward the prolate spheroid axis in 
the equator, then goes near this axis towards the poles, and the returns from the poles toward the equa- 
tor moving beneath the surface of the droplet. 


The possibility of division of active liquid droplets, due to Marangoni effects produced by chemical 
reactions in the surface of the droplet was studied also. 

Mathematical models of the coupling between surface chemical reactions and viscous flow dynamics 
allowed to study the instabilities produced by effects due to oscillating enzymatic reactions (Smith- 
Sorensen, 1980) (Smith-Sorensen and Castillo, 1980). 


Recently, active droplet instabilities, motion and division was studied with the intention of construct- 
ing mathematical models and physical models of prebiotic protocells (Whitfield and Hawkins, 2016) 
(Zwicker et al, 2017) (Weber et al, 2018). 


In 1932, Bugenberg de Jong (Shaw, 1992) found that dilute solutions of polymers coalesced into drop- 
lets when shaken. In these sets of droplets, known as coacervates, organic substances the were in solu- 
tion before the operation, are highly concentrated (Evans and Wennerstrém, 1999). 

The droplets they form can concentrate in its interior certain solutes, as biological cells can do, for 
example with potassium ions. 
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However, these colloidal droplets do not have neither an equivalent to the plasma membranes of the 
cells nor a sustained network of chemical reactions that could be considered as a precursor to metabo- 
lism. 

Later came the discovery of proteinoids and their study (Fox, 1973, 1976 and 1978). Although mainly 
passive from a chemical kinetic point of view, and able to grow by accretion from a solution, divide 
and show other processes suggesting some of the processes seen in cells. 

From a thermodynamic point of view, these phase separations are possible when polymer molecules 
lower their energy enough when they come closer together and overcome the entropic effects related 
with mixing with the solvent. 


In a minimal model of two main components proposed by physicists and biologists from Dresden 
(Germany) collaborating with people from other places (Zwicker et al, 2017) the droplet material con- 
stitutes a phase made of component B that separates from the solvent and forms a droplet (Figure 8.4, 
to the left in green). 


Figure 8.4 (Adapted from Zwicker et al, 2017) 


This droplet structural component B spontaneously degrades to the second component, A , following 
an irreversible reaction (represented by the arrow (1) in Figure 8.4), and diffuses out of the droplet (the 
environment is represented to the right in Figure 8.4, in blue). 

The substance B is formed again by a reaction A+C — B+C"’ (represented by the arrow (2) in Fig- 
ure 8.4) that involves a fuel C that supplies energy and whose concentration remains constant in the 
environment. 


The level of detail of the mathematical modelling of these processes is higher than the one we em- 
ployed in the present article, because it is a distributed parameters approach, while our approach was a 
lumped parameters one (with the only exception of the Marangoni effects that were included in the 
relative elongation rate (7.36)). 


Chemical reactions inside the droplet in the interface and its surroundings are coupled with mass dif- 
fusion and surface tension effects and included in a nonlinear mathematical model of the correspond- 
ing anisocoric open physical-chemical system. 

Suitable partial differential equations of reaction-diffusion in the bulk of the droplet and in the envi- 
ronment are coupled with equations describing surface phenomena at the droplet’s boundary. 
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Combining analytical technics to study stability and bifurcations with numerical simulations of the 
dynamics, they could predict active droplet grow and division, following the stages of this last process 
(Figure 8.5) (Zwicker et al, 2017). 

So, simple chemically active droplets can grow to the size of contemporary cells and after that, spon- 
taneously divide. 


Figure 8.5 (Adapted from Zwicker et al, 2017) 


Some people believe that this mechanism might have enabled liquid droplets to evolve into living cells 
in early Earth’s primordial soup. 

Others think that membranes appeared first to stabilize the incipient and labile self-duplicating auto- 
catalytic systems (Eigen, 1971) already present in a primeval soup before the appearance of living 
cells, through a relative isolation from the brutal changes in the environment. 


A mechanism was proposed (Goldacre, 1958) to explain the spontaneous formation of double layers of 
lipo-protein membranes from single layer films in the interface between air and aqueous solutions. 
Tubes would be formed first, whose diameters are of the same numerical order of the diameter of 
many biological cells, and then a partition in little bags could happen due to mechanical perturbations 
from the environment. ° 

But this leads us to the last topic of the present work: the emergence of prebiotic systems and its pos- 
sible relevance to solve the problem of the origin of life. 


(8.4) Prebiotic systems and the problem of the origin of life. 


The title of Schrédinger book “What is life?” could be restated, from an evolutionary perspective and 
considering the need to fully understand the enormous set of data available in molecular biology, as 
"what are cells and where do they come from?" 

A complete enumeration of the known properties of the different cell types is possible, but it is 
doomed to be insufficient as new data is added. 

Nevertheless, some aspects seem to stand out: life as we know it is based on cells, these cells have 
DNA-based genome that gives rise to a variety of enzymatic or structural proteins, these proteins con- 
nect through networks of chemical reactions with other molecules or ions, and signal from this metab- 
olism together feedback the genome. Information coded in the genome is translated to the protein 
world, and through the cell reproduction the genomic information is transmitted to the offspring. The 


8 The action of the wind drags the lipo-protein monolayers, located at the interface between bodies of 
water and the atmosphere. When the drag is strong enough, tubes like those proposed by Goldacre could 
be formed. This could happen in our time. However, the very forces that gave rise to the tubes will destroy 
them, unless a proper stabilizing mechanism could be implemented. 
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coded information suffers certain changes, mutations, and cell line that interacts with a possibly 
changing environment participates in Darwinian evolution. Cells incorporate high energy and low 
entropy nutrients to anabolic networks that generate structures and send low energy and high entropy 
excreta from catabolic networks that partially destroy the results of anabolism, while anabolic net- 
works are continuously rebuilding these structures in an endless effort until the cell dies or divides. 

To what extent are these and other known properties of cells necessary to an understanding of life on 
Earth? Which of these properties are related with life’s origin in a prebiotic world? If there is a rela- 
tion, in which temporal order they appeared during the passage from the inanimate to the animate 
world? ° 

In part, all this is matter of speculation, due to the uncertainties in reconstructing past prebiotic scenar- 
ios and given a definition of life itself. 

Some people think that to try to identify a minimum set of functional moduli under laboratory condi- 
tions could be a more rewarding approach to try to answer the posed question in relation with the cell. 
Two methods have been used to try to do this. 

One method begins with an already existent cell, establishing tentative criteria to eliminate progres- 
sively parts in the genome while retaining what are considered its main functions. 

The other method consists in establishing procedures to assemble a cell from what are tentatively con- 
sidered essential functional moduli. 

Applying this second method of cell synthesis it could be conceived that something like prebiotic sys- 
tems or protocells would be found along the way. 


Contemporary research on the dynamics of open systems capable of growth, with the synthesis in the 
laboratory of a prebiotic system as the ultimate goal, requires the consideration of open physical- 
chemical systems capable of spontaneously originating spatial patterns of concentration and possibly 
asymmetric osmotic flows, that could modify the balance of mechanical forces and at the right mo- 
ment lead to the division of the pre-biotic system (Szostak, Bartel, and Luisi, 2001) (Hanczyc and J. 
Szostak, 2004) (Macia and Solé, 2007 a) (Macia and Solé, 2007 b) (Rocheleau et al, 2007) (Zhu and 
Szosta, 2009) (Schrum, Zhu and Szostak, 2010) (Whitfield and Hawkins, 2016) (Zwicker et al, 2017) 
(Seyboldt and Jiilicher, 2018). 


The kind homogeneity hypothesis, made in sections (4), (5) and (6) of this paper, considerably simpli- 
fies the analysis but leaves out of the scope of this analysis the questions related to the genesis of pos- 
sible mechanical instabilities that lead to a spontaneous partition of the system. 

As consequence, the necessary relationship between the surface and the volume of the system must be 
entered independently. 

However, in section (7) the study of mechanical instabilities that could lead to a spontaneous partition 
of the system was introduced, albeit limited to determining bifurcation points beyond which the spher- 
ical shape of an active droplet becomes unstable. 

Only in the previous subsection (8.3) of this discussion references were made to more complete math- 
ematical models of shape instability, leading to the splitting of chemically active droplets. 


° For example, a world of protein polymers first: primitive networks of catalyzed chemical reactions ap- 
peared first in certain locations in the oceans, until the first container emerged, and then in their relatively 
protected interior a precursor of metabolism developed including high molecular weight polypeptides, 
until a molecular code emerged and allowed a qualitative new degree of physical-chemical stability and 
individuality to the prebiotic system, opening the door for a Darwinian evolution. Or a world of RNA first, 
with the advantages given by its autocatalytic properties, as many people think (Powner, Gerland and 
Sutherland, 2009) (Nelson and Cox, 2017). 
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People who consider that some type of active droplets could have been the precursors of life on Earth, 
think that these droplets may have evolved to constitute prebiotic systems themselves first, and true 
living cells later. 


So, let us consider the conditions to have a prebiotic system. 
As already said in section (5), it is generally admitted now that three subsystems must interact to con- 
stitute a prebiotic system: 


(1) A network of chemical reactions, precursor of metabolism. 


(2) A macromolecular element, a precursor to the genome, that carries information. 
It catalyzes its own production and that of the structural elements of the system - the container 
- regulating the network of chemical reactions. 


(3) A container, which delimits the prebiotic system of its environment. Able to grow, it regu- 
lates the transport of mass to and from the environment, ensuring necessary conditions for the 
operation of the chemical network and, where appropriate, the replication of the information- 
carrying subsystem. 


-In principle, the initial container could be a coacervate droplet made with low molecular weight pol- 
ymers that could be part of the hypothetical primordial soup, forming membrane-free containers able 
to concentrate chemical substances and thereby foster chemical reactions. Alternatively, the container 
could be some type of membrane bounded vesicle compatible with the composition of this soup, and 
able to concentrate chemical substances and foster chemical reactions. 

These containers could have been produced by some of the mechanisms seen at the end of subsection 
(8.3) above, or other mechanisms compatible with known laws of physics and chemistry and with the 


past geochemical environment. !° 


-The components of the primordial soup, inorganic and organic molecules and ions could be produced 
in thermal vents in the deep of oceans or due to ultraviolet radiation, electric discharges in the atmos- 
phere, in relation with volcanic activities or even in connection with nuclear geysers (Wachtershau- 
ser, 1990) (Morowitz, 2008) (Nelson and Cox, 2017) (Maruyama et al, 2018 and 2019). 

Cell biochemistry is characterized by macromolecular catalyzers, most of them proteins. But these 
proteins, to be synthesized from the corresponding monomers require other high molecular weight 
macromolecular catalysts. And is reasonable to assume that these were not available in the primordial 
soup. 

A way out of this chicken-and-egg like problem is provided by the properties of catalysts formed by 
transition metals in the periodic table of elements with low weight ligands. 

These metals could be present in the early oceans and coastal regions in several oxidation states ligat- 
ed with simple substances and gave rise to more complex molecules driven by thermal or chemical 
gradients. Some of these molecules are catalysts. Laboratory evidence shows that these compounds 


10 Qld rock structures, the stromatolites, are built by thin layers originated by microorganisms and were 
dated approximately as 3.5 billion years old. However, the geochemical conditions in the Earth, let us say, 
4.0 billion years ago, in relation with the primordial soup are difficult to establish. In turn that the traits 
of the geochemical environment must have influenced both the composition of the primordial soup and 
the chemical reactions that could have occurred there. For example, the assumptions about the composi- 
tion of the early atmosphere have changed quite radically from the fifties (when the Miller and Urey ex- 
periments were designed and carried out) to the present time. 
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can catalyze a variety of organic chemical reactions, including the formation of monomers and their 
anabolic polymerization reactions, that would have been necessary during the early stages of biogene- 
sis. These initial networks of chemical reactions could originate autocatalytic loops like a primitive 
version of the tricarboxylic acid cycle (Morowitz, Srinivasani and Smith, 2010). 


-And then we have the third essential component, the macromolecular element precursor to the ge- 
nome, able to catalyze its own production and that of the structural elements of the system (the con- 
tainer), regulating the network of chemical reactions and receiving control signals from an operating 
machinery, already of a proto-biochemical nature. 


Now, let us consider the consequences of each of these three subsystems (container, macromolecular 
element precursor to the genome, and network of chemical reactions, precursor of metabolism) on the 
transport processes and chemical reactions involving chemical species already existent in the abiotic 
world on Earth. 


From one side, there is a physical image of the world that integrates, with a minimal set of fundamen- 
tal laws, processes from the scale of elementary particles to the scale of stars and galaxies, using as 
few constituents as results possible. Sometimes this perspective enhances a reductionist view of the 
cosmos, like that implicit in a possible misunderstanding of the poetic phrase of Hubert Reeves (1994) 
"We are nothing more than stardust, each element that makes us contains carbon, oxygen and other 
elements that have all been, one day, synthesized in a star, which, dying in a thousand fires, spread the 
necessary seed for our existence”. From at least a practical point of view, a full reduction of living 
systems to physics is not feasible, even if it could be done in principle. In any case this is something 
that is very debatable from the epistemological point of view (Anderson,1972). 


From another side, a somewhat inverse approach that, not ignoring the fundamental laws of physics, 
begins by recognizing that those same microscopic laws can give rise to different structures. The or- 
ganization of a physical-chemical or biological system can change abruptly when a parameter crosses 
a threshold, and systems of different natures may behave in an analogous way. 

Different collective motion emerges from disordered microscopic evolution, when only differences in 
boundary conditions modify the effects due to fluxes of energy and matter. Considering all this as well 
as other arguments (Medawar, 1974 and 1978) (Bunge, 1979 and 2003) (Mayr, 2004), an image of the 
world as a manifold of emergent systems, integrated hierarchically is obtained. 

A higher emergent level can be thought of as a boundary condition within which the lower emergent 
entities are integrated and which is in turn constrained by the levels above. 

In this sense, a genome or one of its precursors can be considered as a generator of boundary condi- 
tions that constraint physical-chemical morphogenesis either in a cell or in a protocell. 

There is always an identifiable meaningful part in a genome who has this function. It is given physi- 
cally by a certain linear sequence of nitrogenated bases, whose precise ordering cannot be explained 
by physical forces alone, because from the energetic point of view the different sequences do not dif- 
fer significantly. 

As stressed by Polanyi in a classical article published in Science (Polanyi, 1968), “it is this physical 
indeterminacy of the sequence that produces the improbability of occurrence of any particular se- 
quence and thereby enables it to have a meaning- a meaning that has a mathematically determinate 
information content equal to the numerical improbability of the arrangement.” 

But information measured this way is not the actual information content in the genome, which is in 
fact the prerequisite to initiate and to control the development of mechanism that will produce con- 
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straints to the basic physical-chemical processes of life. This information cannot be derived from the 
laws of physics or chemistry, it pertains to another, emergent, level of complexity. 

A different matter, at least in principle, is how the genetic code could have originated. 

Some people, thinking in the framework of evolutionary biology, assumes that the origin of the code is 
a sequence of random chance variations and selection steps acting on some kind of precursor that 
appeared in some form once, or many times. An estimate of the probability of such events is difficult 
if it is at all possible. Depending on the method of calculation, it could end in an estimated probability 
well below Borel’s limit. |! 

The experimental approach through a fast and faithful reconstruction of the spontaneous appearance of 
a proto-code and its progressive modifications by trial and error, mimicking the effects of natural se- 
lection, seems to be out of the question: for now, impossible. 


Other people think, as the Spanish surrealist artist Salvador Dali: "And now, the announcement of 
Watson and Crick on the DNA. This is for me the real proof of the existence of God." 

This quote appears at the beginning of the essay by Crick “Of Molecules and Men” (Crick, 1966), as 
an antithesis to what the author tries to make clear in this book: how irrelevant he considers religious 
beliefs to a deep understanding of the universe. 

However, the uncertainties that surround the problem of the origin of life open a wide spectrum of 
views and beliefs. 

In any case it seems advisable to remember that ultimately, the fidelity of nature (the assumed invari- 
ance of its laws without which science as it is now would not be possible) is nothing more than a met- 
aphysical postulate. As Max Born once said (Born, 1949): “Yet there are metaphysical problems, 
which cannot be disposed of by declaring them meaningless, or by calling them with other names, like 
epistemology. For, as I have repeatedly said, they are “beyond physics” indeed and demand an act of 
faith. We must accept this fact to be honest...Faith, imagination and intuition are decisive factors in 
the progress of science as in other human activity.” 


Whatever it is, if the morphogenesis in living systems, done by physical-chemical processes, occurs 
under the control of a coded symbolic description, then as Howard Pattee says (Pattee, 2015) the dis- 
tinction between matter and symbol extends from the very origin of life, throughout biological evolu- 
tion and beyond. 
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(9) Appendix 1: A brief review of some topics in dynamical systems theory 
(9.1) Deterministic dynamic systems 


(9.1.1) Generalities 
Considered as a deterministic dynamical systems, the dynamic equations of an open reaction- 
diffusion system with volume dependent of its composition, that include partial differential equations, 
describe the movement of a point x (the state of the system), '* in an abstract space X of infinite di- 
mensions, the phase space of the system. This movement follows an evolution equation that can be 
written like this (Attle-Jackson, 1989) (Logan, 2008): 

dx» 
— = F(x;c) 
dt 
Here F represents a nonlinear operator that may include partial derivatives relative to the Euclidean 


spatial coordinates, c = (as Crysis a) is a vector whose components are fixed parameters and f rep- 


resents the time instant. In the so-called reduced order models of reaction-diffusion systems, the 
space X is of finite dimension and the evolution equation is a set of nonlinear ordinary differential 
equations. When time doesn’t appear explicitly in F , as is the case considered here, the evolution 
equations are called autonomous. 

Each set of possible values of the parameters c = ( Cys C55 5G; ) can be considered as a point belong- 
ing to a k-dimension Euclidean space C=R*: c= (cs ChiscriCp ) eR‘ Wecall C= R* the param- 
eters space. 

Given enough regularity, for each initial condition x, there is a unique trajectory x(t ; aye) . The set 
of phase points that belongs to a trajectory determine an oriented curve in phase space name orbit. 
The family of all these curves is the phase portrait of the dynamical system for the chosen set 
C,,Cy,+..C, Of parameter values. 

When a change in some parameter or parameters (called control parameters) produces a qualitative 
(topological) change in the phase portrait, it is said that a bifurcation has occurred. 

For example, a stable equilibrium (fixed) point x (c) is destabilized and a stable limit cycle of oscilla- 
tion appears (a kind of local bifurcation of the dynamic type) or new equilibrium (fixed) points appear 


at a certain distance from the already existent equilibrium points in phase space without changing the 
type of pre-existing fixed points (a kind of global bifurcation of the static type). 


The bifurcation diagram of a deterministic dynamical system is the partitioning of the parameters 
space into regions where the phase portraits are qualitatively (topologically) equivalent. 

A set included in phase space, such that if a point of a trajectory belongs to the set, then all the points 
of the trajectory belong to this set, is called an invariant set. (A fixed point, the points of a cyclic tra- 
jectory and in general the points belonging to a set of trajectories is an invariant set). 

Sometimes invariant sets, regular enough to be considered as a refinement and a generalization of the 


regular curves and surface of ordinary Euclidean space EF, can be found in the phase space of a dy- 


namic system. This kind of sets, known as finite dimension invariant manifolds locally resemble an 


2 The state of this distributed parameters system is composed of several concentration fields defined in 
spatial regions that in general change with time if the system grows or if the system self-consumes. 
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Euclidean space EF’, (the same E,, everywhere, son is the dimension of the manifold), but globally it 


might be more complicated. (For example, a torus and a Klein bottle are two-dimension manifolds). 
An attractor is an invariant set such that all trajectories passing through nearby points tend to that set. 
(Stable fixed points and stable limit cycles are attractors). 

As will be seen below in relation with the mathematical foundation of the methods of reducing the 
number of state variables and generate the so-called reduced order models, certain higher dimension 
invariant manifolds can be attractors. 


In relation with rest points (fixed points or steady state points) their stability or instability are of inter- 
est in the case of either growing or self-consuming open chemical-physical systems. 

This stability can be approached using two methods of Liapunov (Luenberger, 1979) (Attle-Jackson, 
1989) (Beltrami, 1997) (Jordan and Smith, 2007) (Logan, 2008): the so called first method that implies 
a linearization in a neighborhood of the rest point, and the second method that implies the use of Lia- 
punov functions. 

Now, we mention some results related with autonomous systems of ordinary differential equations 
(ODE) because the mathematical models of open reaction-diffusion systems constructed in the present 
article, even those including bulk diffusion, are of this type. 


Consider an open set GC R” and a reduced order model given by an autonomous system of ODE 


x = f (x) being f (x) a vector field with continuous partial derivatives in G . 


Let Q2G by a bounded open set such that x, is the only rest point of the system of ODE there. 


Suppose that we find a smooth function V(X):Q5R' that has a minimum only at x, and 


“v( x(r)) < (along every solution curve. This function is called a Liapunov function in Q for the 
t 


autonomous system of ODE (Beltrami, 1997). 


If there is a Liapunov function the rest point is locally stable in the small, in the sense that the orbits 
that begin in a suitable small neighborhood of x, remain in Q. 


If “v(x(1) <(Qin -{x,} the trajectories that begin in in a suitable small neighborhood of x, 
t 
tend to X, (asymptotic stability of x, in the small) which is now an attractor. 


If “v(x(1) > Qin some neighborhood of the fixed point x, (excluding x, because “V(x) =0) 
t t 


the fixed point is unstable (it is a repellor). 


The following invariant set theorem (Luenberger, 1979) (Beltrami, 1997) if suitably adapted, is use- 
ful to study the stability of the steady states of the anisocoric model in subsection (4.2.2) “Scale effects 
in the anisocoric generalization of Bertalanffy’s model of open system”. 

a = oe 7a : 
We begin with an autonomous differential system = =F (x) with the vector field f (x) having 
t 


continuous partial derivatives in an open setG CR” . 
- Let V (X)be a scalar function defined in G with continuous first order derivatives relative to the 
components of x. 


Let Q, ={¥eG/V(x)<s} be a bounded region and assume that if xeQ 


s? 
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then VV (x ) of (x) <0 The symbol V represents the gradient operator and « indicates a scalar prod- 


_— 
IL 
a 


uct of the vector fields, so VV (x) -f (x) = $27) He (x) 


Furthermore, let S = {x eQ /VV (x) -f (x) = 0} and let I, be the largest invariant set included in 


the set S . 
-Given the above-mentioned restrictions, every orbit beginning in a point of ©, tends to I, as time 


increases. 


If the largest invariant set I, has only one element, the rest point x, of the open system, all orbits 
beginning in Q, (as well as all the trajectories corresponding to geometrically equivalent systems) 


tends to this rest point. If I, is a cycle in phase space, all orbits (and trajectories) tend to this attracting 


limit cycle as time increases. 
In relation with the qualitative behaviour of the orbits, two observations are useful for the case of an- 
isocoric models with scale effect. 


~ P 

An invertible change of variables y = Ax (or in components y, = a ij (i SL 2ae P)) transforms 
j=l 
dx oe fe aN s dy a ee : (= 7 7( 4-1 =. > < rem | P 
rr = f (%¢) into the system Fe UH) being h(y) =Af(A y3é), being A the inverse 
t t 

matrix of A= ae ‘ 
Both systems are equivalent in the sense of having orbits with the same pattern of topological (geo- 
metric) properties (Attle-Jackson, 1989). This result is of use in subsection (4.2.2) to ease the con- 
struction of a Liapunov function. 


(9.1.2) Nonlinear finite dimension positive systems as mathematical models of anisocoric open 
systems 

Now, let us restrict our considerations to the special case of the open physical-chemical systems with a 
finite vector 7 of state variables that motivated the present research. 


Let us define R? such that if n € R? , for every nj ( pa leks p) which is a component of 71, then 


n, 20. The closed set R? is the only subset of R’ with physical meaning when n, ( ; ee p) 


are the number of moles of the components of a reaction-diffusion system. 
The mathematical models of finite dimension reaction-diffusion systems belong to the class of posi- 
tive dynamical systems in the sense that the nonnegativity of the components of the state vector are 


preserved during its motion in phase space. Due to this, the largest possible phase space is a R? . 
In the mathematical models of open anisocoric reaction-diffusion systems, described with lumped 


parameters, whose volumes depends of their composition and has scale effects, the nonlinear vector 


eek i a Ta 
field f (7i,c) of the corresponding system of ordinary differential equations " a (7,¢) that de- 
t 


scribes their dynamics, is not defined for all 7 € R’? because it has a singularity in the points of R’ 


where the system volume is zero, either being unbounded in their neighbourhood or not having a 
limit in these singular points. 
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This is always the case of the origin 7 =0. However, as we can see in the case of the generalization 
of the Bertalanffy model, anisocoric and with effect of scale, the volume of the system is zero for other 


points of the space of states: it is so for every point (n,n, ) = (n,,0) with n, 20. 


If ( R? ), is the open set formed by the interior points of R? , then the field f (i, C ) that describe the 


dynamics of open physical-chemical systems is smooth in (R? ) and the theorem of existence and 
“ ae | brine gs : 
unicity of the solutions of ae eg (vi, é) can be applied to any initial point 1 (0) € (R i" ) : 
t v 
However, in case of a self-consuming system, from any initial condition the corresponding future tra- 
jectory verifies lim, ,,,. 7 (t) =0 and 0 is not in the domain of f (i 3C ) ; 
If the focus is not in kinetic properties of the trajectories properly, but in the geometric patterns of the 


orbits of the system and in their stability properties, it is possible to make modifications in the vector 
field without altering the orbits (Nemytskii and Stepanov, 1960) (Roxin, 1972). 


Instead of working with f (a,é ) it is possible to work with the bounded field 


1 


nF@ey -f (a, c) where ||f ( 7i,c) 


| represents, as usual, anorm. The new field has the same 


direction and the same zeroes as the original vector field. So, it generates the same orbits. 
This result is a regularization theorem (Roxin, 1972). 
Instead of leaving outside the domain of the system possible limit points like the origin of a self- 


consuming system, the bounded field can be modified so that the new domain includes the whole R? . 
This result is a consequence of Vinograd’s theorem (Roxin, 1972). 


a°(ii,(R?)"] 
+a" (i(k?) } 


a point 7 € R” and the closed set (R ‘a )" which is the complement of (R i ) relative to R”. 


We define the scalar function g (i ) = 


me 


where d (i.(RP )'} is the distance between 


Then g (7i) is positive, bounded between 0 and 1, different from zero in every point of (R . 


) and zero 
in the points belonging to (R . \ : 


Then, the vector differential equation — 


dt aaa ay} ° f(a,2)= f(a 7i,C) has the same 


geometric patterns of orbits in (R i ) that the original vector differential equation (a local dynamical 
L 


d Peat 2 
system relative to R’ ) a =) (7i,c) but has been extended regularly to the whole R’ (Nemytskii 
t 
and Stepanov, 1960) (Roxin, 1972). 
All the points of (R’ ) are rest points of the new dynamical system. 


The regularization and Vinograd’s theorems can be applied to the geometrically equivalent system that 


results after an invertible change of variables y = Ai. This is done in section (4.2.2) “Scale effects in 
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the anisocoric generalization of von Bertalanffy’s model of open system”. 


(9.1.3) Linear finite dimension positive systems as mathematical models of anisocoric open sys- 
tems 

When there is not a scale effect and the chemical reactions are, from a kinetic point of view, chains 
of reactions of the first order, as in the biophase of Perret and Levey or in our model of biophase as an 
extension of Bertalanffy model of open system, the evolution equation is now defined in the whole 


dni x .y., 
R? and can be written: a = K(é)n 
t 


Here K (é ) is a kinetic matrix of a class known as Metzler matrices, characterized by having 
nonnegative off diagonal elements and real eigen value A(é ) such that the real parts of the other 


eigenvalues are strictly smaller than /, (é ) , So this last eigenvalue is dominant (Luenberger, 1979) 
(Attle-Jackson, 1991). 
It has a corresponding real eigen vector V, (é ) that is a dominant eigen vector in the sense that, as- 


ymptotically for t + +00, we have: “i(t;¢) ~ exp [4, (€): t| ‘Vy (E) 


The expression exp[ A (é ) : t| Vv, tc ) , where A, (é ) is an eigen value and v, (é ) is the correspond- 
ing eigenvector of the kinetic matrix, is called a normal mode of evolution. 

The state vector n (t5¢ ) can be decomposed in a linear combination of normal modes of evolution. 
Almost always the eigen values of K (¢ ) are all different (no degenerate eigen values), either real or 


in pairs of complex conjugates ones (Luenberger, 1979) (Hirsch, Smale and Devaney, 2004) and then 
the corresponding set of eigen vectors is a linearly independent set. If this is not the case for a certain 


matrix K F (é ) , introducing one norm in the space of square matrices '’ we find that there are matrices, 
with eigen values all different, as near as we want from the matrix K . (é ) : 


The measurements obtained in experiments always have a certain uncertainty, so the components of 
any matrix related with measurements will be uncertain and at least in principle it is possible to dis- 


card cases like K : (¢ ) (Attle-Jackson, 1991). However, there are certain situations in chemical kinet- 


ics that is more convenient to describe with matrices that have degenerated eigen values (Starzak, 
1989). 


The eigenvalues of K (e ) are continuous functions of the parameters. The first eigenvalue that even- 


tually changes its sign when a control parameter crosses a critical value is the dominant eigenvalue. 


Before the crossing of the mentioned threshold value of the parameter, n (36 ) tends to the origin of 
phase space R?. After crossing the threshold, ‘i(t 3 ) tends to become parallel to Vv, (é ) growing 


asymptotically as an exponential with the now positive coefficient A (é ) F 


(9.2) Random dynamic systems 


Although the theory of dynamic systems and one of its branches, bifurcation theory, ignore fluctua- 


13 All the norms are equivalent in finite dimensional Banach spaces like this space of square matrices 
(Griffel, 1983) (Hirsch, Smale and Devaney, 2004). 
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tions, as Hermann Haken (Haken, 2003) (Haken and Portugali, 2016) has repeatedly emphasized fluc- 
tuations are crucial precisely at and in the close neighborhood of the bifurcation thresholds. 
If random fluctuations in the parameters are incorporated in the mathematical model, a random dy- 
namical system is obtained, described by a random differential equation: 

ae F (%c+ d¢(t)) 
dt 


Here 0¢ (t ) represents the random fluctuations of the system parameters, X(t ) is a stochastic process 


’ OKs feat 
that represents the random state of the system at each instant of time and oe is its stochastic deriva- 
t 


tive (van Kampen, 2001) (Gardiner, 2004) (Scott, 2013). 
In the limit of large numbers of particles involved in random events in a representative volume ele- 
ment of an open reaction-diffusion system, the random state can be expressed as the sum of a deter- 


ministic component x(t ) and a random component OX(t ) é X(t) = x(t) + 0X (t) 


This approach is known as the linear noise approximation (van Kampen, 2001) (Scott, 2013). 


In the chemical kinetics of open systems considered as prebiotic systems, we assume that the term 


x(t ) is a not chaotic solution of the usual full scale deterministic equations of the reaction-diffusion 
system, and the fluctuating term Ox(t) is described by a probability density distribution 
II (Ox, t) centered upon x (t) : 

It is possible to imagine that IT (Ox, t) moves along the deterministic term x(t ) : 

When the phase space is of finite dimensionn , the fluctuation term Ox(t) is composed by n sca- 
lar random processes 6X, (t), 6X, (t),..., 0%, (¢). 

The corresponding multivariate probability density T1( Ox, (4); OX, (@); Se Oae (t) :t) is in general 


Gaussian and is obtained solving a Fokker-Planck equation that describes the evolution of the proba- 
bility density of the fluctuations. (Gardiner, 2004) (Scott, 2013). 


Away from the bifurcation points the width of the probability distributions of the state variables are in 
general small relative to their mean values. This allows the estimation of the deterministic values of 
the parameters of the system, from measured data, by well-known statistical methods (Gershen- 
feld,1999). 

When the linear noise approximation can’t be applied, there are more powerful methods that in princi- 
ple allow us to distinguish between chaos, almost periodicity and random noise, as well as to identify 
possible attractors of the dynamic system. 

These methods are based on embedding techniques and the construction of information dimensions 
(Gershenfeld,1999). 


As stressed by several authors (van Kampen, 2001) (Gardiner, 2004) (Scott, 2013) random fluctua- 
tions can modify the bifurcation scenario. 

However, under conditions of applicability of the linear noise approximation, a deterministic analysis 
can explain the stability behavior. 

This applies to the mathematical models of open reaction-diffusion systems considered in this article. 
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(9.3) Methods of reducing the number of state variables 


(9.3.1) Inertial manifolds 

When the main interest is in asymptotic dynamics, we find that often it is possible to substitute a 
mathematical model of an open reaction-diffusion system based in nonlinear partial differential equa- 
tions by a mathematical model of the same system based in nonlinear ordinary differential equations 
that show a similar behavior. The origin of this similarity in the asymptotic behavior of dissipative 
partial differential equation (DPDE) and certain ordinary differential equations (ODE) is due to the 
following fact: the asymptotic evolution of the solutions of these DPDE occurs near and exponentially 
approaching to a finite dimensional invariant attracting manifold in phase space, the so-called inertial 
manifold. For each trajectory that begins in a neighborhood of the inertial manifold, there is a trajecto- 
ry included in the inertial manifold such that the distance between these two trajectories tends to zero 
exponentially. These and other properties of the inertial manifolds are described by several authors 
(Constantin et al,1989) (Malinietski, 2005). 

Now, this inertial manifold corresponds to or is very well approximated by the solutions of a suitable 
set of nonlinear ordinary differential equations: with the inertial manifold theory, it is shown that cer- 
tain dissipative PDEs have the same asymptotic behavior of a finite dimensional ODEs system, known 
as the inertial form corresponding to the considered inertial manifold. 

Furthermore, the exponential approach of the solutions of the DNPD to the corresponding inertial 
manifold (solutions of the ODE) is in many cases fast enough to allow the dynamic study to be done 
entirely on the inertial manifold, that is, using ODE and reducing computational effort. As conse- 
quence, if an inertial manifold exists, in principle at least, a reduction of the dimension of the state 
vector of a complex system can be done, preserving the main properties of the original model and 
attaining a small approximation error in the framework of the simplified mathematical model. 

Unlike other asymptotic methods, in principle at least, inertial manifold theory could be applied to 
simplify the analysis of a dynamic system at all points of its parameter space: a reduction by inertial 
manifold. 

However, the construction of the exact inertial forms that correspond to a given system of DNPD is 
often not feasible. So, in practice only approximate systems of ODE (approximate inertial forms) can 
be constructed to reduce the dimension of the phase space of a complex system. 


Besides inertial manifold theory (a global method), there are two general procedures to reduce a high 
dimension dynamic to a low dimension dynamic: the center manifold theory (a local method) and the 
slow manifold theory (a global method) (Nicolis, 1995). 
In both cases a few dominant mode amplitudes can often be identified, so that the behavior of the sys- 
tem is essentially determined by these dominant modes. 


(9.3.2) Centre manifolds, deterministic and stochastic 
In the center manifold theory, the stability of the equilibrium solutions (fixed points) x (c) in phase 
space X is studied in parameters space C = R‘ varying systematically some parameter values (con- 


trol parameters). 
The system is expressed in each neighborhood of an equilibrium solution as the sum of two terms, one 


linear J (x(c))-(x-x(c)) and the other non-linear h(%(c);x-x(c)). This last one vanishes 
faster (at least quadratically) than the linear term when the state of the system x tends to the equilibri- 


um point x (c) in phase space. So, we have: 
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dx 


4 = 4(¥(6))-(x-x(6)) +n(()sx-¥(6)) 
For a nonlinear system of ordinary differential equations of order n, X = R”. 


The corresponding fixed points x (c) are called hyperbolic (or non-degenerate) when all the eigen- 


values of their Jacobian matrix J (x (c)) have non-zero real parts. This occurs almost always in pa- 


rameters space. 
According to the Hartman-Grobman theorem, also called principle of linearized stability (Attle- 
Jackson,1989) (Nicolis, 1995), the behavior of the dynamic system in a neighborhood of a hyperbolic 


point X(c) is topologically equivalent to the behavior of the linearized approximation 


& J (3(c))-(x-3(0)). 


dt 


However, for certain points c,- (bifurcation points or critical points) in parameters space C=R* the 
corresponding equilibrium point x (ci) in phase space X = R” can be non-hyperbolic (or degener- 


ate): they have some eigenvalues with zero real part. In case of complex conjugate eigenvalues, there 
is at least a pair of complex conjugate eigenvalues laying on the imaginary axis of the complex plane. 
So, the bifurcation points (critical points) in parameters space, where the real part of one or in general 


a few eigenvalues of the linearized system in a neighborhood of x (c) change its sign, are identified. 
In the local development “ =J (x(c)) : (x —X (c)) + h(x(c) 3x x(c)) of the dynamic equations, 
t 


the nonlinear term h(x (c);x-x (c)) can‘t be left aside when one seeks to build a dynamic system 


topologically equivalent but easier to analyze than the original. 


In a neighborhood of these critical values of the parameters, where the equilibrium points lose their 
stability, the change of stability is often preceded by the critical slowing down of a few mode ampli- 
tudes (the dominant ones) that slave the others in such a way that (perhaps after a short transient) the 
evolution of the state of the system occurs in a low dimension manifold (which exactly at the critical 
point is the so called center manifold in bifurcation theory) (Nicolis, 1995) (Haken, 2003). 

The other modes of evolution are given as functions of the dominant ones, called order parameters in 
Synergetics (Haken, 2003). 

This reduction in the description of the dynamic system to a dynamic in terms of the dominant modes 
of evolution is known as center manifold reduction. It generates a reduced order model. 

By near identity transformations of state variables (Nicolis, 1995) (Kuznetsov, 1998), in a neighbor- 
hood of the bifurcation point it is possible (although not always easy) to build a dynamic system that 
gives a local description of the dynamics topologically equivalent to the original system, but easier to 
study. The nonlinear differential equations of these topologically equivalent dynamic systems are 
known as normal forms and allow a classification of the local bifurcations. 

In principle the normal forms related with center manifold reductions are asymptotically exact, unlike 
the inertial forms related with the reductions to inertial manifolds accessible in practice, which are 
asymptotically approximate. 

The concept of a normal form as an asymptotic local description of the bifurcation dynamics of a sys- 
tem in the neighborhood of an equilibrium point, is not exclusively related to the classical approach 
initiated by Poincare, the method of multiple time scales or the more recent variant developed by 
Arneodo, Coullet and others (Nicolis, 1995). Normal forms can be obtained also by the averaging 
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methods (Suarez-Antola, 2919). 

Local bifurcation theory has several advantages to cope with stability problems in reaction-diffusion 
systems: many stability issues can be understood as a change in the stability condition of a steady state 
operating point, including the birth of limit cycles of oscillation (Hopf bifurcations) either stable (su- 
percritical) or unstable (subcritical). 

But other stability problems seem to be outside the span of local bifurcation theory. For example, the 
derivation of the boundaries of instability around a stable steady state (operating point) in case of hav- 
ing at the same time at least another stable steady state with the possibility of an undesired run away 
from the operating point to this other steady state. 

When random fluctuations can't be neglected, it is possible to resort to a random center manifold theo- 
ry, which allows a stochastic approach to bifurcation theory. 

Random normal forms appear, in general as nonlinear Langevin equations like this one, in terms of the 


random distance a (t) from a steady state and a random phasew (t ) : 


<ai(t)=«-a(t)-7-@ (0) +N, (0) <H(t)=0() +N, (1) 


Here N, (t) and N, (t) are random noises, «is areal and y is real and positive. 


Without noise, the corresponding deterministic normal form describes a supercritical Hopf bifurcation 
that occurs when x changes its sign from negative to positive. In presence of the noise this formalism 
describes a stochastic supercritical Hopf bifurcation. 

If the linear noise approximation is valid, usual statistical methods can be applied to estimate the nu- 
merical values of parameters, in this last case «and y, from noisy data. 


In a neighborhood of the bifurcation point, when linear noise approximation usually fails, but N ‘ (t) 


is a Gaussian white noise of zero mean and known variance o”, it is possible to derive a formula 
for the probability density function (pdf) of a (1) averaged over on cycle of oscillation (Suarez-Antola, 


2019). 


, : kK : : 
The only parameters that appear in this formula are —; and a (unless numerical factors), so if o7 
Oo Oo 


is known, the other parameters could be estimated from averaged random data. Besides, it is possible 
to apply the more powerful methods framed in embedding techniques and information measures. 


(9.3.3) Slow manifolds 

The slow manifold theory, when applicable, allows us to identify a few dominant degrees of freedom 
in phase space, even when the equilibrium solutions are far from bifurcation points in parame- 
ters space. A so-called slow manifold, of low dimension, must be identified in phase space. Further- 
more, this slow manifold must be stable in the sense that any trajectory located in a neighborhood of 
the manifold, after a short transient, should approach to and remain near this manifold. It is called slow 
because after this approach to the manifold, the state of the system changes slowly in comparison with 
the initial transient. 

According to this definition, the center manifold defined in , on which the relevant part of the dynam- 
ics lies, is an example of slow manifold, because fast modes of evolution are slaved by slow modes in 
a neighborhood of a bifurcation point. But not every slow manifold is a center manifold (Nicolis, 
1995): it is not necessary to have a bifurcation point of the dynamic system in the set of parameters 
values considered in slow manifold analysis. 


Moreover, a slow manifold is not necessarily formed by orbits of the dynamic system (a slow mani- 
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fold is not necessarily an invariant manifold like the inertial manifold is), although some orbits must 
pass close enough to the slow manifold for a certain interval of time at least (Tihonov et al., 1970). 


Now we are going to summarize some relevant aspects of Tihonov’s approach to singular perturbation 
theory, which gives a foundation to identify and characterize slow manifolds. 

Let us suppose that the equations of evolution of a deterministic dynamic system can be written as 
follows (Tihonov et al., 1970) (Lin and Segel, 1988) (Malinietski, 2005): 


dz dy 
an ry 3 ay : 
os (z, y3c) 7 f (z yse) 


The state variable is given now by the pair of vectors (2; y) interconnected through the preceding 
evolution equations. As before, c represents a point in parameters space and the new positive parame- 
teré is small enough and in the limit will equal zero. This means that z(t sc) varies much faster 


than y (t; c) when Z is not too close to F (z, y3 c) =(). The initial conditions of the given trajectory in 


phase space are z(t3c) = Z, and y(t3c) = Vis 
Now we consider the degenerate dynamic system whose evolution equations are obtained from the 


original ones putting ¢=0: 


d 
F(z, y;c)=0 at (exe) 


In general, the solution of the equation F (z,.93c) =0is a set of branches 
z=h,( y3c) (k =1,2,...,m) that do not intersect and can be classified as stable branches 


z=h, ( y; c) and unstable branches z = h, ( y3c) A branch h( y; c) is stable if all the eigenvalues of 
OF 
the operator a (h ( y3 c) V5 c) have negative real parts (Tihonov et al., 1970) (Malinietski, 2005). 
zg 


These stable branches are the stable slow manifolds mentioned above: given initial conditions belong- 


ing to the region of attraction of the stable branch, the trajectories Z (t ; c) and y (t ; c) first approach 


to the stable branch and then remain near it. 


Since z(t3c) changes much faster than y(t;c) as they approach the stable branch, it is possible to 


introduce the following simplification of the mathematical model, known as the inner approxima- 
tion, keeping the slow variable fixed at its initial value: 


d 
e— = F(% 95) 


This simplified equation is solved with the initial condition z(t3c) =Z, and an inner solution 


Rivet (t; c) can be constructed. 


Once near the stable branch, a second simplification of the mathematical model is obtained, known as 


the outer approximation, substituting z =/, (y3c) in the equation of the slow variable: 


—= f(h,(y3¢), 3c) 
This new simplified equation is solved with an initial condition y, that verifies z) =h, (%: e); and an 


outer solutions Y,,,., (t ‘ c) and Z,1er (¢; c) =h, ( Lone (t; c);c) 
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The inner and outer solutions can be combined to obtain a valid approximation including the fast re- 
laxation from the initial conditions to a neighborhood of the stable slow manifold, and the slow 
movement near the manifold (Tihonov, 1970; Lin and Segel, 1988). 

The slow manifold theory summarized above may be applied to study transients, including nonlinear 
oscillations, because in a complex system there is a hierarchy of widely separated time scales. 

When the variables of interest (variables of reference) belong to a certain time scale, and there is an 
attracting slow manifold, it is possible to simplify the dynamic analysis applying the following two 
principles to link scales of different orders of magnitude: 


1-The variables belonging to processes with time scales at least an order of magnitude greater than the 
reference time scale, can be considered as frozen. This corresponds to setting the values of the slow 
variables to build an (approximate) inner solution (in Tihonov’s sense) for the fast variables. 

2- The variables (fast) belonging to processes that evolve with time scales at least an order of magni- 
tude smaller than the variables of reference (slow), after a short transient (produced in the so-called 
inner time scale) can be considered as relaxed to equilibrium with the reference variables (evolving in 
the so-called outer time scale). This corresponds to the building of an (approximate) outer solution 
(in Tihonov’s sense) for the slow variables. 


Using these principles, both the slow manifold and a few dominant degrees of freedom of evolution, 
each with its corresponding amplitude (a function of time), can usually be identified. The dominant 
degrees of freedom evolve more slowly and tending to slave the evolution of the other degrees of free- 
dom. So, an infinite dimension dynamic or at least a finite dimension dynamic with a big number of 
dimensions, is reduced to a low dimensional one: a reduced order model is obtained. 

Often slow manifold theory can be applied far from bifurcation conditions, to further reduce the di- 
mensionality. 


(10) Appendix 2: The scale effect in increasing objects of constant shape 


To define objects of constant shape, let us define first the meaning of similar regions as composed by 
the union of star-shaped regions. 


A star shaped region is a closed region R, of R® such that there exists an interior point O that can be 
connected to each point of R, by a segment of straight line included in R, . 
In spherical coordinates centered in O, a star-shaped region can be described by 


O<p<f (9, 2) 0<@<a O<@<2z where the function f (9, 2) is continuous and piecewise 


continuously differentiable (Lin and Segel, 1988). 


Two star-shaped regions RO and Re? are called similar if, after possible rotations and translations 
they can be described by the equations 0<p< f (0,9) and O<p< Ff (0,9) with 
f? (0,9) =L- fo (0,9) for a positive constant L. 


Between the surfaces and the volumes of two similar star-shaped regions there are the relations 


[RO ]=2-v[R] and s[R? |=2-s[ Rr, ]. 


1 (8 | RY | 
=—-| —=—~= | and relative to a fixed region R, 
Livy | RS ) 


(1) 


it can be seen that 


As consequence —=—7S 
v| (2) 
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RO 
the quotient Fool decreases when the length L increases. 
R, 
7 - 2 
s|RO} (vi rR}? S| RO a 
Eliminating L we obtain: 7 = Fat or S$ Re = AVE. : (v [rR |p 
s|RO} LvjpR® (7\3 
ee te") 
2 
So, for similar star-shaped regions we obtain: S [R, | =U: (Vv [R, 13 


These results remain valid for more general regions that can be decomposed in a finite number of star- 
shaped regions that share at most part of their border areas, when the same length L is applied to ex- 
pand the component star-shaped regions and to produce an object of the same shape. 

Since the volume of an object of constant shape increases as the cube of a characteristic linear dimen- 
sion, while the surface increases with the square of this linear dimension, the quotient between surface 
and volume decreases with the increase of this linear dimension. In the case of biological organisms 
considered as open physical-chemical systems, as the volume increases the surface tends to be too 
small to supply the interior with nutrients and to excrete waste products. 


(11) Appendix 3: Reciprocity theorems of continuum mechanics 


Let us consider a linear elastic body in equilibrium under a system of volume forces with force densi- 


ties (forces per unit volume) F and surface forces per unit area (stresses) T. The corresponding 
displacement field produced in the body by the volume forces and surface stresses is S . 


The region filled by the body in R’ is represented by B and its boundary surface by OB. 
The following reciprocity theorem of Betti and Rayleigh (Segel, 2007) relates two systems of forces 


and displacements (F, He ; 5) and (F, T, 155 ) , applied on the same elastic body: 


[se °5,)aV +[(0 5, ds 7 [, (#2 °5JaVv +[,(% 5, )ds (A3.1) 


The symbol « indicates a scalar product of the vector fields. 
In the derivation of this equality two theorems are employed. One relating the stress and strain tensors 
with each system of applied volume and surface forces: 


[, (Aes )av+]. (Tos, )as = (> Oe: eV (A3.2 a) 
i, j=l 


[, (Bes av +. (Tes,)as = bs Oj, 4 }a¥ (A3.2 b) 
i, j=l 


In this last equalities o, ,, and ¢ are the components of the cartesian stress and strain tensor fields 


Lij Lij 
for the first system of applied forces, while o,;,; and €,,, are the corresponding components of these 


tensor fields for the second system of applied forces. 
The other theorem is another reciprocity theorem that relates the same tensor fields: 


3 3 
i 2.6 ef 28, JOY (A3.3) 


i j= i.j= 
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If we choose o, .. =1 (i =J= z) and the other components of this stress tensor equal to zero, from 
(A3.3) and (A3.2 a) it follows: [42.0 =] (Fes, )av+[_ (705, jas (A3.4) 


If E is Young’s modulus and v is Poisson’s modulus of a linear and homogeneous elastic solid, the 


displacement vector field 5, is given by (Love, 1927): 

a iC a nu = 

8, (uyz)=a[ee v(x +y-2)} (A3.5) 
From (A3.4) and (A3.5) and representing by V (B ) the volume of the body, we obtain for the aver- 


&€,..:dV inthe z direction: 


1 


age strain (e, a = 


(Sec) Bay qy al -v-(x-F,+y-F,))avV+[ (2-7, -v-(x-T,,+y-T,,))-ds} 


(A3.6) 
The reciprocity theorems (A.3.2) and (A.3.3) can be extended to steady varying deformations in gen- 
eral linear hereditary materials. These generalized versions of the reciprocity theorems involve rates of 
deformations and rates of displacements and accelerations are neglected (Rabotnov, 1980). 
If stresses are related through linear relations with rates of strain with symmetric tensorial coefficients 
(Rabotnov, 1980), it is possible to write the following formulas that are the viscous analogous to the 
elastic formulas related with Betti’s and Rayleigh’s reciprocity theorem: 


-~ O._ ~~ OO. 3 O 

j{F “Ss, Jv, [7 +, Jas -j(a, So, av (A3.8 a) 
fe. MO ae FOr 3 0 

j{F “Ss lave (7, “Saas -)(Seu Sa, Jv (A3.8 b) 


i, j=l 
- 0 3 0 
(> Or, Sa, av = (> oe Say, jv (A3.9) 
i, j=l i 
If we choose again o, .. =1 (i =j= z) and the other components of this stress tensor equal to zero, 


from (A3.8 a) and (A3.9) it follows: 


0 = 0. = O. 
|, ae dV = AG “25, )av +f[F +, Jas (A3.10) 
If a steady state of viscous plastic flow is established after a yield point has been reached in the mate- 


rial, it is possible to describe the deformation giving to Poisson modulus the value v = 5 (incom- 


pressibility) and substituting Young’s modulus E by 3-77 where 77 is the corresponding viscosity of 


the flow (Goodier, 1936). 

Following Goodier’s suggestion, Nadai (Nadai, 1937) introduce a relation between the stresses and 
strain rates for viscous creep in an isotropic and homogeneous elastoplastic solid when elastic defor- 
mations can be neglected. 

In this case the material behaves as incompressible and a formula for the velocity of displacement can 
be derived, analogously to the derivation by Love of formula (A3.5): 


86 


August 2020 


One 1 . 1 S = 
— S,(X,¥,Z) =) 7'e, - = Xe. tye, A3.11 
Fa, (s9.2)=504 oa v-2,)} (A3.11) 
= , O 1 0 
Defining the average rate of strain eg ears } 6 ay 3 from (A3.10) and (A3.11): 
ar /~ V(B) /8ar 
Oe ete | zF =4-(x-F, +y-F );-av+] z-T, Ste +y-T,,) |-d5 
Ot 1,zz 3-7-V(B) B 1,z 4 ix y ly aB 1,z 5 ix y Ly 


(A3.12) 
Youg (Young, 1939 a) proposed to use this last formula to describe elongations in a mathematical 
model of prolate spheroidal (ellipsoids of revolution) cells framed in the continuum mechanics of vis- 
cous bodies. 


Besides (Young, 1939 a), considered the possibility of having the density of volume forces F given 
by the gradient of a scalar potential y: F=V y (A3.13) 


If (A3.13) applies, the volume integral that appears in formula (A3.12) can be transformed in a surface 
integral over the boundary of the body (the unitary vector field 1 represent the outward normal in 


each point of the boundary 0B: 


JEFF (Rte k) av =f, [ed 0-5 (0g eiity-d, 07) } yeas (A3.14) 


In the case of the active fluid droplets considered in section 7 above, the combination of a porous pol- 
ymer matrix and a multicomponent solution that fills the interconnected space of pores can be mod- 
elled as an equivalent elastoplastic-viscous and incompressible material. Assuming a quasi-steady 
state viscous flow it should be possible to apply (A3.12), as a convenient approximation to study the 
deformation of an initially spherical active droplet. 
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